{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing a custom kernel in GPyTorch\n",
    "\n",
    "In this notebook we are looking at how to implement a custom kernel in GPyTorch. As an example, we consider the [sinc](https://en.wikipedia.org/wiki/Sinc_function) kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we start, let's set up some training data and convenience functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iter = 2 if smoke_test else 50\n",
    "\n",
    "# Training data is 100 points in [0,1] inclusive regularly spaced\n",
    "train_x = torch.linspace(0, 1, 100)\n",
    "# True function is sin(2*pi*x) with Gaussian noise\n",
    "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * math.sqrt(0.04)\n",
    "\n",
    "# Wrap training, prediction and plotting from the ExactGP-Tutorial into a function, \n",
    "# so that we do not have to repeat the code later on\n",
    "def train(model, likelihood, training_iter=training_iter):\n",
    "    # Use the adam optimizer\n",
    "    optimizer = torch.optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters\n",
    "\n",
    "    # \"Loss\" for GPs - the marginal log likelihood\n",
    "    mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "    for i in range(training_iter):\n",
    "        # Zero gradients from previous iteration\n",
    "        optimizer.zero_grad()\n",
    "        # Output from model\n",
    "        output = model(train_x)\n",
    "        # Calc loss and backprop gradients\n",
    "        loss = -mll(output, train_y)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        \n",
    "\n",
    "def predict(model, likelihood, test_x = torch.linspace(0, 1, 51)):\n",
    "    model.eval()\n",
    "    likelihood.eval()\n",
    "    # Make predictions by feeding model through likelihood\n",
    "    with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "        # Test points are regularly spaced along [0,1]\n",
    "        return likelihood(model(test_x))\n",
    "\n",
    "def plot(observed_pred, test_x=torch.linspace(0, 1, 51)):\n",
    "    with torch.no_grad():\n",
    "        # Initialize plot\n",
    "        f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "        # Get upper and lower confidence bounds\n",
    "        lower, upper = observed_pred.confidence_region()\n",
    "        # Plot training data as black stars\n",
    "        ax.plot(train_x.numpy(), train_y.numpy(), 'k*')\n",
    "        # Plot predictive means as blue line\n",
    "        ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')\n",
    "        # Shade between the lower and upper confidence bounds\n",
    "        ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)\n",
    "        ax.set_ylim([-3, 3])\n",
    "        ax.legend(['Observed Data', 'Mean', 'Confidence'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A first kernel\n",
    "\n",
    "To implement a custom kernel, we derive one from GPyTorch's [kernel class](https://docs.gpytorch.ai/en/latest/kernels.html) and implement the `forward()` method. The base class provides many useful routines. For example, `__call__()` is implemented, so that the kernel may be called directly, without resorting to the `forward()` routine. Among other things, the `Kernel` class provides a method `covar_dist()`, which may be used to calculate the Euclidian distance between point pairs conveniently.\n",
    "\n",
    "The `forward()` method represents the kernel function and should return a `torch.tensor` or a `linear_operator.operators.LinearOperator`, when called on two `torch.tensor`s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FirstSincKernel(gpytorch.kernels.Kernel):\n",
    "    # the sinc kernel is stationary\n",
    "    is_stationary = True\n",
    "\n",
    "    # this is the kernel function\n",
    "    def forward(self, x1, x2, **params):\n",
    "        # calculate the distance between inputs\n",
    "        diff = self.covar_dist(x1, x2, **params)\n",
    "        # prevent divide by 0 errors\n",
    "        diff.where(diff == 0, torch.as_tensor(1e-20))\n",
    "        # return sinc(diff) = sin(diff) / diff\n",
    "        return torch.sin(diff).div(diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now already use this kernel. We therefore define a GP-model, similar to the tutorial on exact GP inference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the simplest form of GP model, exact inference\n",
    "class FirstGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super().__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = FirstSincKernel()\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By using the convenience routines from above, the model can be trained and evaluated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/arthus/miniconda3/envs/torch-toying/lib/python3.8/site-packages/torch/autograd/__init__.py:130: UserWarning: CUDA initialization: Found no NVIDIA driver on your system. Please check that you have an NVIDIA GPU and installed a driver from http://www.nvidia.com/Download/index.aspx (Triggered internally at  /opt/conda/conda-bld/pytorch_1603729096996/work/c10/cuda/CUDAFunctions.cpp:100.)\n",
      "  Variable._execution_engine.run_backward(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAADGCAYAAADWg+V4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvvklEQVR4nO29e3yT9dn4//4kbZOeC6WFQjlWzrSUgggDUVDQIeJE+AniDh4eh45Nf4/scT7uUdC5Oac4N6d7MafORwd4QN3Q+QgbjrMCCsixAkIptIWWntskTfL5/nEnN2matEmTpmn5vF+vvGju3IerJfd1X9f1uQ5CSolCoVAYOlsAhUIRHShloFAoAKUMFAqFC6UMFAoFoJSBQqFwoZSBQqEAwqAMhBBmIcTnQoh9QoiDQogV4RBMoVBEFhFqnoEQQgCJUso6IUQssBW4X0q5MxwCKhSKyBAT6gmkpk3qXG9jXS+VyaRQdDHCEjMQQhiFEHuBc8AGKeVn4TivQqGIHCFbBgBSSgeQL4RIA94TQoyRUh7w3EcIcQ9wD0BiYuL4ESNGhOPSCoUiCPbs2VMupczw9VnIMYMWJxTiMaBeSvmMv30mTJggd+/eHdbrKhSKthFC7JFSTvD1WThWEzJcFgFCiHjgWuBIqOdVKBSRJRxuQhbwFyGEEU25vCWlXB+G8yoUiggSjtWE/cC4MMiiUCg6kbAEEBXdl6amJoqLi7FYLJ0tiiIIzGYz2dnZxMbGBnyMUgaKVikuLiY5OZlBgwah5Zcpoh0pJRUVFRQXFzN48OCAj1O1CYpWsVgspKenK0XQhRBCkJ6eHrQ1p5SBok2UIuh6tOf/TCkDRdRTXFzMTTfdxNChQ8nJyeH+++/HZrMB8Nprr7F06dJOlrAlSUlJPrcbjUby8/MZPXo0Y8eOZeXKlTidzlbPdfLkSf761792hJjNUMpAEXZKSkq46qqrKC0tDflcUkrmzZvHd77zHb7++msKCwupq6vjkUceCYOkvrHb7R127vj4ePbu3cvBgwfZsGEDH330EStWtF7oGyllgJQy4q/x48dLRdfg0KFDQR9z7733SoPBIO+9996Qr79x40Z55ZVXNttWXV0te/bsKevr6+Wrr74q586dK6+77jo5bNgwuXz5cimllHV1dXL27NkyLy9Pjh49Wq5Zs0ZKKeXu3bvltGnTZEFBgZw1a5Y8e/aslFLKq666Sj788MNy2rRpcvny5XLgwIHS4XBIKaWsr6+X2dnZ0mazyWPHjsnrrrtOFhQUyKlTp8rDhw9LKaU8ceKEnDRpkpwwYYL8+c9/LhMTE33+Pt7bjx8/Lnv27CmdTqf85ptv5NSpU+W4cePkuHHj5LZt26SUUl5xxRUyJSVFjh07Vq5cudLvft74+r8Ddks/96VSBopWCUYZmM1miVax2uxlNpvbff3nn39ePvDAAy225+fny3379slXX31V9unTR5aXl8uGhgY5evRouWvXLvnOO+/Iu+++W9+/qqpK2mw2OXnyZHnu3DkppZRr1qyRd9xxh5RSUwaeymvu3LnyX//6l77fXXfdJaWUcsaMGbKwsFBKKeXOnTvl9OnTpZRS3njjjfIvf/mLlFLKF154IWBlIKWUaWlpsrS0VNbX18vGxkYppZSFhYXSfZ9s2rRJ3nDDDfr+/vbzJlhloNwERdg4ceIEt912GwkJCQAkJCSwePFivvnmm3afU0rpMxjmuX3mzJmkp6cTHx/PvHnz2Lp1K7m5uWzcuJGHHnqILVu2kJqaytGjRzlw4AAzZ84kPz+fX/ziFxQXF+vnvPXWW5v9vHbtWgDWrFnDrbfeSl1dHdu3b2fBggXk5+fzwx/+kJKSEgC2bdvGokWLAPjud78b9O8IWk7Hf/zHf5Cbm8uCBQs4dOiQz/0D3S9YVJ6BImxkZWWRkpKCxWLBbDZjsVhISUmhT58+7T7n6NGjeffdd5ttq6mp4fTp0+Tk5LBnz54WykIIwbBhw9izZw8fffQRDz/8MLNmzeLmm29m9OjR7Nixw+e1EhMT9Z/nzp3Lww8/zIULF9izZw8zZsygvr6etLQ09u7d6/P49kTwT5w4gdFoJDMzkxUrVtC7d2/27duH0+nEbDb7POa5554LaL9gUZaBIqyUlZWxZMkSdu7cyZIlS0IOIl5zzTU0NDTw+uuvA+BwOHjwwQf5wQ9+oFsgGzZs4MKFCzQ2NvL+++8zZcoUzp49S0JCArfffjvLli3jiy++YPjw4Zw/f15XBk1NTRw8eNDndZOSkpg4cSL3338/c+bMwWg0kpKSwuDBg3n77bcB7Ym+b98+AKZMmcKaNWsAePPNNwP63c6fP8+SJUtYunQpQgiqq6vJysrCYDDwv//7vzgcDgCSk5Opra3Vj/O3X8j48x868qViBl2H9gQQw01RUZGcM2eOvOyyy+SQIUPk0qVLpcVikVJK+eqrr8oFCxbI2bNnNwsgfvzxxzI3N1eOHTtWTpgwQe7atUtKKeWXX34pr7zySpmXlydHjRolV61aJaXUYgbufdy8/fbbEpCffvqpvu3EiRPyuuuuk3l5eXLkyJFyxYoV+nZ3APFXv/qV35iBwWCQY8eOlaNGjZJ5eXnyN7/5jR6oLCwslLm5ufKKK66QP/vZz/Rz2Gw2OWPGDJmXlydXrlzpdz9vgo0ZhL2fQSCofgZdh8OHDzNy5MjOFkPRDnz933VoPwOFQtE9UMpAoVAAShkoFAoXShkoFApAKQOFQuEiHA1R+wshNgkhDrvGq90fDsEUCkVkCYdlYAcelFKOBCYBPxJCjArDeRUKQMvs80zxtdvtZGRkMGfOnE6UqvsRsjKQUpZIKb9w/VwLHAb6hXpehcJNYmIiBw4coLGxEdAyDvv1U1+xcBPWmIEQYhBap2Q1Xk0RVr797W/z4YcfArB69Wq9KAigvr6eO++8k8svv5xx48bxwQcfAFofgCuvvJKCggIKCgrYvn07AJ9++ilXX3018+fPZ8SIESxevJjOSL6LNsJWqCSESALeBR6QUtb4+FwfrzZgwIBwXVYRQR54APzU6LSb/Hz47W/b3m/hwoU8/vjjzJkzh/3793PnnXeyZcsWAJ588klmzJjBK6+8QlVVFRMnTuTaa68lMzOTDRs2YDab+frrr1m0aBHuzNcvv/ySgwcP0rdvX6ZMmcK2bduYOnVqeH+5LkZYlIFrFPu7wJtSynW+9pFSrgJWgZaOHI7rKi4d8vLyOHnyJKtXr2b27NnNPvvkk0/429/+xjPPaBP9LBYLRUVF9O3bl6VLl7J3716MRiOFhYX6MRMnTiQ7OxuA/Px8Tp48qZRBqCcQWt3mn4HDUsqVoYukiFYCeYJ3JHPnzmXZsmV8+umnVFRU6NullLz77rsMHz682f7Lly/3W+prMpn0n41GY4e2OusqhCNmMAX4LjBDCLHX9Zrd1kEKRbDceeedPProo+Tm5jbbft111/H73/9e9/u//PJLoANLfbsp4RivthXokF7ab+06TUp8DL1TzGSlxpORbMJoUG27L1Wys7O5//6WaSz/8z//wwMPPEBeXh5SSgYNGsT69eu57777uOWWW3j77beZPn16s+YlipZEdQnzqs3Hqbde1OYxBkFmiklXDn1SzaTGBz4+ShE8qoS56xJsCXOXantmd0rOVlk4W2XhS6oASDQZdeWQlWqmd4qZuBiVZa1QBEuXUga+qLc6OHG+nhPn6wEwCEHPpDiyUsz0STWTlWqmZ2KcmgqkULRBl1cG3jilpLzWSnmtla/OVANgijXQR1cOmgVhjjV2sqQKRXTR7ZSBL6xNTk5VNHCqokHf1iMhlqy0ePqkaNZDryQTBhWcVFzCXBLKwBeVDU1UNjRx6KyWLBkXYyAz2aQHJrNSzSSaLtk/j+ISRH3bXdjsToorGymubNS3pcbHkpV60b1QS5uK7owKu7dCdWMTR0pr+fToeVZ/XsSLm47x1q7TbC48z9dltdRZVdZapCgtLWXhwoXk5OQwatQoZs+e3Sy9OFC2bNnC6NGjyc/P58yZM8yfP9/nfldffTWXWgdvZRkEgd0pOVPVyJmqi9ZDsjlGC0qmaa5FZrK5W1sPz20I/gZsjf9/5rA295FScvPNN/P9739fH1Syd+9eysrKGDas7eM9efPNN1m2bBl33HEHAO+8807wQndTlGUQIrUWO4Vltfz76HnWfH6aFzcdY+2uIv7tsh5qLU2dLWKXZ9OmTcTGxrJkyRJ9W35+PlOnTuWnP/0pY8aMITc3V5+N6K9E+eWXX+att97i8ccfZ/HixZw8eZIxY8YA0NjYyMKFC8nLy+PWW2/VeyeAVgg1efJkCgoKWLBgAXV1dQAMGjSIxx57jIKCAnJzczly5AgAdXV13HHHHeTm5pKXl6ePh/N3nmhBWQZhxjMx6gvXNrf10CfVTN+07m89hJsDBw4wfvz4FtvXrVvH3r172bdvH+Xl5Vx++eVMmzYN8F2ifPfdd7N161bmzJnD/PnzOXnypH6ul156iYSEBPbv38/+/fspKCgAoLy8nF/84hds3LiRxMREfv3rX7Ny5UoeffRRAHr16sUXX3zBiy++yDPPPMPLL7/ME088QWpqKl999RUAlZWVbZ4nGlDKIALUWuzUWmopLNPm5bnTqvukxtPXFaBMNqu06mDZunUrixYtwmg00rt3b6666ip27dpFSkpK0CXKmzdv5ic/+QmglUvn5eUBsHPnTg4dOsSUKVMAsNlsTJ48WT9u3rx5AIwfP55167Tq/Y0bN+ruDECPHj1Yv359q+eJBpQy6ATCYT2UlJSwcOFC1q5dG9KU467A6NGjffr2rdXVtKdE2d/o95kzZ7J69epWr+N5DeljjHxb54kGVMwgSnDHHjYXBhZ7eOKJJ9i6dSuPP/54J0kcOWbMmIHVauVPf/qTvm3Xrl306NGDtWvX4nA4OH/+PJs3b2bixIntusa0adP06ckHDhxg//79AEyaNIlt27Zx7NgxABoaGtpcxZg1axYvvPCC/r6ysrJd54k0ShkEQU3FOV548HZqLpwPaHsouK2HL05Vsn5/CS9v+YaXt5wgzmRGCMFLL72E0+nkpZdeQghBfHx82K4dbQgheO+999iwYQM5OTmMHj2a5cuXc9ttt5GXl8fYsWOZMWMGTz/9dLutpHvvvZe6ujry8vJ4+umndaWSkZHBa6+9xqJFi8jLy2PSpEl6oNAfP//5z6msrGTMmDGMHTuWTZs2tes8kaZLlTB3Nu/8bjk7PlzD5BsWMv8ny/1ur6k4x+u//E++98hzpPTMCKsMNRXn+NuqX/PV9o00WS3EmcxMm3UDK558irxhg0gKc9akKmHuunTrEubO4r/m5GG3WfX329evZvv6lr6fe7swGEBKPnnjD82URjhISc/EnJCE3WYlJs5Ek81Kg4zjs1Inn5WeINkcQ980V+xBZU0qgkApgwD4+V82Nnsax5rM5E6ZyfQFd7HprZf17W6k0wlcVA4xcSaeXr8/bPLUVlXwrTmLmDT7VnZ+tLaZe1JrsXO0tJajpdrKRaxRkJls1pOislLjVc2FwifqW9EKnua+59PYbrNiTkiiX87IFtsBDAYjTqdDVxpz73koqGu15Vrc8djF4NQtP36s1X2bHC2zJlNcNRdZquZC4UG4WqW/AswBzkkpx4TjnNHAJ2++yDcHdvPJG3/w+zSuraoAIZq5EU6nFudoslowJyS1uLl93fjua61/+RkulJ3h5vse4b0Xn+yYuENjEzWNTQFbD76WyhTRTXtigWEJIAohpgF1wOuBKINoDyB6xwjc+DP3vYN6wmBk+PgpJKWlY22oa/Ykh+YBx88/WefzWqBF0b2Dle0l2KBmSnwsfVPNZBlr6ZPeg8yMdAwGtfjUFZBSUlFRQW1tLYMHD272WYcHEKWUm12j1boF3jGCmDgT5vhEfvjrV5vt58+NcDTZ6Nm7X4ub2FcgEkAIA1I6W8ghpQw67uDvpve0cgJRLm7r4bhwMrKulFNnSzAKgdFw8WVQ1kLUYjab9QzMQIlYzKArjVfzjtjbbVbqbFZ2rF/T7EYKxI3wxF8gUhgM7Nn4gU+lEBNnIm/qLKYvuIsXHry9xU3uffN7uxtFR/Zhb7Lp+werXJqkgf3VJp+fea5cXAoVm92dsOUZuCyD9d3BTQB4dcVSDuz4p74yEAiB3GDvPP8YOz5aizE2DkeTjck3LKS2spyUnhnUVl1g/5aP9X3dS5STb1gI0GqOA0L4llUICq6+oYUCmnvPQ2GPRXjWXLgDlKrmIrpozU1QyqAVvGMB/pYU/d1gvkz2V1csJaVnRjMLwh1TcH9WdvoEdZUVJPVI5/j+z4NSSK3hdmHCFYcIhGRzjG459EmNp3eyiRijij10Fl0y6ShM3/+AcN+03hF8X+6CvyVFX6sGvvz01pYFvQONbtkCVUip6b0pP3uqmbsRazKTkJzG0HGTuWreD/y6MB2FVrFZx9dlWu2+0SDolWTS28n1TY0nNUFZD9FAuJYWVwNXA72EEMXAY1LKP4dyznvvhQ8+7k/2UAv9h1kYMNxC3xwrcabwp0+7b9o3nlrGuaLjur/9vUeea3VJ0V+MwF/GYnuSj4JRSE6ngyk33qa7G0IYsNusjJ40XVdGbeUldDQOp6SsxkJZjQVOa9sS4oz0STXTJ8VM37R4MlNMmGJUK/tIE7W1Ca+8Ar97pZ5vDpuouaDpLINRkjXYqiuHASMs9Blgw9DO742/JUQ3wmDg2Y8Pt9je1jKdv6d5e/10T9di87rXOPjZJv5r1Xre/f0Kny5Ha65IV0AISE+M02MPfVLNpKtBOGEhIjGDYAg2ZlBVHsPpoyaKjpo5fdRMUaEZS72mAeLMTrKHXlQOA4Zb6JFpp7XvjadbsOmtl9m/bUOrSsH7ie6vYMkTX4HCWYvva6ZE2lPQFMi1uyNxMQZ9xkUf1yshLmq93KilyysDb5xOKD8bS9ERM6cLzRQdNXPmmAl7kxaYSkqza8phuIWBIyz0H24hIfliEMLzhkJKdny01mc03vuJHkwykq+nc3KPXs1u5GBu7GAToS4FPFvZ90k1k5GkgpNt0e2UgS/sTVDyjWY9uF/niuKQUjMRMrJtlJ99B+ncAXwG7AO09XdhMJCTN5G6ygoa6qqpqTiHEAZANrtR22v+/9cNuc3W+v3h68b2tmIisUTYVTEaBBnJJj3+kJVqJi0hrrPFiiq65GpCsMTEQv9hVvoPszLlRm3GYmO9geJCE6eOmCk6Gk9j3S3UVd3mOsKKEPvp1a+MqTcNY8SEWHr1beK1x1s+0d34C+a1dTPmXzWb3Rvf1wuYYuJMJPfoRW1lOXabtdWCJndwc8eHa9t17UsJh1NSWm2htPpiBWl8nJE+Kdp0breSiI9TwUlfdBtl4Mv/jk90MnRcI0PHNQKVSAl//fXz7PlXGcIwGemcQEXJtbz3BzMACckOBgx/j6S0RmouWLj++6NITGnuOgSSaejG27R3FzDZbVbiTPE4mmx+b2y/qcsGA/c//1bElwi7Ko02B9+U1/NNeb2+LS0hVlMQrvwH5V5odBs3oTX/29PUfumhHzDqiun6mnt1RQXXf++PnDoST9ERM6eOmCk75eFe9LMxYISFgSMbGTjCQtZgKzEBLou3VsBUuGcruVNm+o34h3tFQuEfd+5Dn1QTvVO0qs0eCbHdcvWiW7sJgazpe+YRNNZWE2cy0y9nBDNvu5fXf/mfJKWdYfLsDCbP1twLS4PQApNHzJw6Ek/hlwns+WcKADGxTrKHWhnoUhADRvhfvfB2K/wVMPla+w/WJenIVmvdnWa5D2jfgbgYg+ZapJh1JdHdU6u7vDLwV/wz956HWiiKslNaZ9q22pOZEyRD8xsZmn/Rvag6H6PFHo6YOXU4nm3rU/n3uh4AJPe0a8phRCMDR2pJUqZ4zeIKxq3wJphjg61KVLSOze7k9IUGTl9o0LclmWK02gtX/KF3ihlzbPeJP3QLN8Hfmv4rK5aSkp7JkV2bW80jcBPMMp3DDmdPmHTr4dRhM+fPaJFrYZBkDbIycKRFe41oJCO7iY5oB9DWkqOyGDoWd/whM8VM7xQTmclm4mKiN/7QrdwEX19uzyfo5nWvsXfzP7BZGzl9dD+Z/XNwNNkQBgPS6dT/7dV3INUVZS2siUAxxnisXszVTMv6GkMz5fDlp8ns+DANgPgkhxZ78Ig/eOY+tJfWLCNo22JQyiI0qhqaqGrQpnUDGISgZ2IsmS4Xo3eKucu0letyyqCt4p/d//wA6XSye8P7AJQVHdM/6zNwKEk90undfwiHd/077Mt0iSlORk5sYOREzbR0OuF8cRynDps5eVhzLzb8tSfSqX0xMrNtDBzVqCuIPoNsGIO0Ov3FFn7xvWsCqo9oj3vhr7BLKRZwSkl5nY3yOhuHztYAWoAyPSmO3smacuidaqJXoglDlCmILuMmtGUOt1ZnUDDjxhZR+M7K33cHJ0+5lMOpw2bqqjWdHGd20n+4hUEjG13uhYXkHm27Sb5+l1uWPtrqaoS/RChjbBwDR4xt9YZ2r9xkDsjhXNHxdmVUXurEGAS9kk26a9E7Rau/6GgF0S0yENtaavP+HJo3B4nWL6eUcKE01mU5aArizHETTof2pUjPsmmWwyhNSfQdYsUYoD3nK5bi/jv89emHmiVCuf+eBqORPRs/8Pk3a6uwy5tLOVW6PcQatSVOt2vREQqiS8cM2mpX7n56eZrL7nr+vKmzSErtGdXJOUJAelYT6VlNjJ+h+Z02q6D4axMnD8Vz6oiZY/sS+GKTa2kzzkn/YRYGuYOTIxtJTfdtPfhajfC2CDw7OX/xr7/r291uhTE2jt98qI0Wd8cnvAu7jLFxWqemADIqFf5pckhKqi2UeGRQuhVEpsuCyEwxkZ7YMTGIqFcG/voMugOFs27/kc9Aom4ud3L9fnuIM0mGjLEwZIz2pXAvbZ48ZKboSDwnD5vZ/H4ajre1qHWPzCYtKDlSUxL9LrMQE+u7kYq3ReBOhDpX/A0XSk432y6dDsZdPVs/h1vhegdknfamNjMqFe2juYLQAtUxBkF6konMZE1JyPpKfnLPD0KeyB21bkJ8fDwWi6XFdveT6lL3T+02QfFxk+5anDxkpuq8lhSjJ0aNbNQtiF/eMSooE98bt8nvqzVb7/5D+GrbhlYzKhWh4y9A+87vlrPzo7X88Ic/5MUXX2z1HF0yZrB3716unH4NloZ6zfR3PancTyRvPM3ZS5Wq8hiXctAUxOnCi2XdKT2tGIy7qLnwd5yOzSD2MWLC5cTFJ3Lq0Jc01FZpsRZhILlHOnVVF5BeU6HUk75zcT8Ax19zk8/O127MZjONjY0tT0AXjRmsWrWKuqoL+nvp8m29FYEvc/ZSJa2XnbQr6xh7pdZv0N4EZ46bdQVx6LPhOB1TtZ2lleKviyiY3gtL/9cp/HIlxtgynHYbCUmp1FWWtzD51dJh5+AduN298X3tByEomD5HD5onJCRw880388wzz7TrOuGaqHQ98DxgBF6WUj7V2v6tWQb+3AM3no0/faEi2P55dcVSzAlDyRryffZsPEfl+YHYLKOx21xNYVIbMSV8RV3VRwwrgBm3TmD731/R26x98sYfLmnXrLNwr5R9sWm9333MZjM2m61NV6FD3QQhhBEoBGYCxcAuYJGU8pC/Y1pTBiUlJSxbtoz333+fhoYG/eavKCnSl8h69O5HTu7l1FZVUPjFNpwOZc62F3da9clD8Xpy1IVSLa3aGCOJTzpOXdV6YIfrdVo/VineyPHO84+x/cM1fjtf/+npR/nTn/5ESUkJ69at83uejnYTJgLHpJQnXBdbA9wE+FUGrZGVlUVKSgoWi4VYj66/3qsEi376FO88/xjS6VQR7BDwTKu+8jvatpoLRh6//SEc9gnUVU1GG4T1gOuIMwjDLvoOucCs2ydgs4oO6VitaE5tVUWLztdNVgsy0cGcux5k7Nix/OEPfwjpGuGwDOYD10sp73a9/y5whZRyqdd+nuPVxp865dvMB5g3bx5ZWVn0HP9tNr2/2m9kuqt3AY5mmidx2UGMAzkRIaYg5UQgB9A6Vve7zMrAEY0MciVG9ejdekNaRfvx/M6/+dSDlJ46xrfmLGLr394MqP9CR7sJC4DrvJTBRCnlj/0d01UmKl3qeGYw2m1Weg+8jNt/9iw7P1pLRWkT35rznO5enD5qxmbVYg/uku5BrrqL/sMsxJmV9RAu/GWCtraK4Kaj3YRioL/H+2zgbBjOq+hkfCVx9csZ4ZHIVc+YyVo7MYdDa0h76pBWtXnykJkD25MAMBgkfYdYGTjKolsQ6VlN1F5QqxO+8F618X7vq1J1wS3zePbZZ0O6bjiUwS5gqBBiMHAGWAjc1vohiq5Aa6PgvDEaIfsyK9mXXSzprqsyasuaroYwuzeksO1vaQAkpdqJNddSWTaTt3+7idsfvlVvCNOdaG051t9n3pWk3u99VaqmpKSElH0I4VtanA38Fm1p8RUp5ZOt7a/chEsTpwNKTsbx3NLf4HRMACYBI12fOkAcZPLsAR3eECaStJYpG2gSkTcxcSZGXj6tWbwsM6ax1VUEN10yAxHg4wMlnL7QSJ3VHgGpFJGieXAyHmPsVDL6LSIx5XrOHE/Rp2XFJzsYONyi110MGB6ehjCRoLWSe8B3argQFFx9g27+t9ZS39vKeODaoSEHEKM2AxHg+jFZANRZ7ZTVWDhXY+Vcrda4UlkMXQdvc7i5mStxNK1nyJgk5v9kPE7nec4Xx3HykNmV9xDP0TfS9W7VvQdY9V4PA0c20mdg+2dtdiStdqCS0ncSkZT6NncD3UgWgEW1MnCTZIohKSOJnIwkfVutpYmyGivnaiyU1Vooq7HSaFMKIhpx+7zu6dY33/cIe7d8zPhrv9NiTLzBAL0H2Og9wMYV12udgiz1Bk4XmvTA5MGdiXz+f6nAxYYwA0dcbAoTSEOYjqat7tbmBO277C+JyP13+WrbhnY31A2WLqEMfJFsjiXZHMtlmRcVRHVjk6YcaqyU1lg4V2vB2tQ1zMqujL/ot7cP7M6pf/be7yCE0FvWtxWcNDcbhqOVdFeUxF4MTh6J59N3euB09ASgZ+8mvY39wJEWsnM0CyTStNbd2lcSkd1mZfSk6XpsYeZt91Jy8mtmLr6PlJ4ZHV6OH9Uxg1CRUlLV0ESpqye+281ocnS/qHVn4h0k8wyMOR32VnPqITxpzTar4Mwxk9ZK7ojmYrhLuo2xTvoN0dwLd1Pa9KymqEiMai1xLpgy/XDEDLq1MvCF0ykpr7dyrsaqzeWrsVBRZ8PZCX+Hrk4wbdA8zWE3MXEm8qbO6rB6kuoKo64cio40T4xKTHXPutAUxIDhFuKTosOKbM/E7W4fQOwIDAahtY9KNjOmn+Z3NjmcnK/VXIsyl4KoamjqZEmjH+8gma/ot9sHtlka2b/lY/1YYTDgaLIFFBBrb+l0arqDvKl15E3VSrodDig9Gcepw/EUHdVcjEOfXXQzM/tbdeXgHqUXaL/JcNJW+/uO4pJTBr6INRromxZP37R4fZulyUFZjUW3HtQKRkt8jY/zjn67feBXVyxlyo23teiQFEhALFzTooxG6Jdjo1+OjW/NuTipu8g9KeuImcOfJ7Jrg/aQiIlz0n+olQEjtLTqAcMtEam7aO+071BRysAP5lgjA9MTGZieqG+rsTRR5upHV1pj4XytFZs9OkzLzsI7SOYv+h1MAZm/AKS/2Q+hEJ/oZPj4BoaP12ZdaN2qYyg6ctF62Pb3NP79rqvnQ5pdVwwd6V6EMpavvVxyMYNw4nRKKuptuvWgxR+sqPBDaHgHIDt7ErXD7qq78LAgzp026Z9nZmuTuvsP15RDvyHtW70IpZOUihl0MgaDICPZREayiVw009Jmd2ruhcvFKKuxUGtRGZSB4Le9F3RqzwpjDGQPtZI91MqUGy+6F6cLTRQd0QKUR79IYPdGrZ09oom+QywMHu1gwPBGBgy3YIo/wxtPtX6jd/bwXGUZRIBaS5NuPZRUWzhXo5Y3fXHm2CH++PBdWBvrWwQgPZOTorFnhZRQXR7D2pUfcHSPnZT02VgbRmJt1NwLY0wDDvsO+uXUMfWmYez46Kfc8djDpKZntGv1wBu1tNhFcS9vlrriD2U1Fi7U2y559+Kd3y1n+/rVwMV03K7Sb9H3DW0ARqA1A7vc9e9YQMt/iDVVMmxcLJn9Kzh99M+cPPJn7Laz7XKHlJvQRfFc3szL1ra5Vy9Kqi+uYFwq6dW+biStPb6B2sryTpIqMNx+/v3Pr2XTWy83W2Y1xyey+L9/yq6P3+Wr7Q+6xv6ZgHzgcpqsl3Nw5+Uc3DkSeNz1OkGT9XMqSuycKx5AnNmCOSEyTwmlDKIEX6sXVQ02XTmUVFsor7PicHY/86G1dfVob3ri9vN3fLi2xXJgnc3KV1s+0bdrDX2tGAy7cTo/03/PWbf/N2ufXYfTOQ5T/FWcOnwtpw734sVlIIQkc4CN/kMt9B9upf8wC/1yrMR2QHq1UgZRTFpCHGkJcYzM0gJTTQ4n52qtlFY36kqiOwQnQ11X9xWFD3XGQ1vHe1szbvcG8LldGAw88Lu39b6Fnr9nZnZPfvzc3a4jtASp2spqTn9t5vRRE6cLzRzdk8jujVqQ2mCUZA2ykj3MqimJYRZs08B0cYGjXShl0IWINRrolxZPP4/kKHdw0q0cymos2KPEegjmhvS1rh7o8Z5R+FmL7+P1X/4nPfv08xmZb885fcUs/Fkz0xfc1cxd8LZyMrIHk5M3sc38geQeDkZNrGfURK2tnDtAWXTUzOlCE8Vfm/lqaxKf/UNTEFv/Ivn881b/xG2iAojdDIdTcr7WSkl1o64kqhs7J7U61HmYbR0fTG2EOzLf3nMaYmIZNDK/WWXmsz+aR11leYuR956NZDsyCOpOkDr9tZkbx2Yxf75aTVC0Qb3V7hF7aORcB2dOhrpUFujxzTsm+Z/CVTDjRvZv+T+fLcXcN/nN9z3Cey8+yczF9/LGr5bpy5vu8X0Z2YMoP3OqWWXm9vWrm3WMdi97dkYL/05fWnS1SV+O1shuopQyoDtcKYPOxemUlNdZ9VHfpdWNVIaxMMv7Jg02IOjv+OkL7uK9F5/0+XQ2xMTiaLLpI+VB89ORksk3LGTW4vuandP7Js8ckMO5ouNkDsih7NSxoH9nY2wcA0eM7bROz+FQBqG2mzwAzAM2h3geRQQxGASZKWbG9k/j+jF9+MGUwSy5Koeb8vtyxeCeDOiZQFxM+78awQYEayrO8cKDt+v+s7/jd3y4VvfjQfPray+cJ3NADg/87m36DLwMp9OBEJrseVNn8a05i6itLNfP6bYg3IN8zxefREpJ2alj+r+tYYyNo0fvfnovw1iTmYIZNzLu6tnNZOuKhKs78qfAMmUZdB+klJTX2XTXoqTaQmVD4IlRwZjKvvx4z+NX/ujmFtO3fSGEgW/NWej3mtrg2SRqqyo4umdrq+eMiTOR1qsPFSVFIATS6dQtjcz+OZw7fVwfLuPv+EjOoex0N8HjAp/ShjIIZryaIjqxNDkorbZw1iM4GUrsob2xgWC6BvtaPXAH+Nw3uTfum75H736MvHxai7Lrr7ZtIHfKTCbNvpXN616j8MsdNNRWdWqOREQyEIUQGwFf0xkekVJ+0ObVXUgpVwGrQLMMAj1OET2YY40M6pXIoF5aYpSUWtVmSZVmPZQGmVYdaBOPQPom+HNDfC0Rupcx3Td5bVU5yWm99H89ey346jvouc09AHjHR2u7/ADgNpWBlPLaSAii6HoIIeiVZKJXkoncbG292209lFRfVBD+mtIGE1vwzEPYvO419mz6u8/uym58JQWFuxeCL9ki1XugI1BJR4qw4st6uFBv01cuSqobm1kPgd5Inr5/rMmMdDha7a4cydZhwYyhi2ZCUgZCiJuB3wMZwIdCiL1SyuvCIpmiWyCEID3JRHqSSe856Wk9DPr9K5RUa9ZDWzdSME/7zmod1pUJSRlIKd8D3guTLIpLhLash9LqRip8xB6Cfdp3hPkeas1DNKPcBEWn48t6sNodzWouSqotEOTTviPM987uRtSRKGWgiEpMMc1LuqWUVDY08eFzdcy59ftcMfv/4x9vvxGxYF0kA5KdhapNUHRZrHYHZdVWfdWipNp/Q5hwlDRHc88F1elIcUljijEyID2BAekJ+rZKV+yhtEbLmiyv1aZlhWreXwoBSaUMFN2KHolx9EiMY1RfrSFMfHw8FsvFisZQzPvukk/gD+UmKLo1JSUlLFu2jPfff5+GhgYSEhK44cab+Ml/P47dlEJptYVztV2/nZxyExSKNsjKyiIlJQWLxYLZbMZisdCrZxpT8y7T97G72sl59nzoDu3kgkUpA0W3p6ysjCVLlnDPPfewatUqSkpKmn0e42PWZp3VrveavFRmXSg3QaEIAO9ZF6VBlnR3NMpNUCgihL9ZF/qczW4w60IpA4WinXinVUPXnnWhlIFCEUa8Z124g5O69dCJ3arbQikDhaID8RWcbLDZdcVQ6prY7a/nQyRRykChiDAJcTEMyUhiSEYScLHuoqS6kbIaC6XV1k5xL5QyUCg6GSEEPRPj6JkYx+i+WtWmp3tRFqFhOEoZKBRRiC/3otHm0GMPZTXhX71QykCh6CLExxkZ3CuRwR6rF9UNTZTWWALKMWgLpQwUii5MakIsqQmxYTlXSBOVhBC/EUIcEULsF0K8J4RIC4tUCoUi4oQ6Xm0DMEZKmQcUAg+HLpJCoegMQlIGUspPpJTu8q6dQHboIikUis4gVMvAkzuBf4TxfAqFIoKEZbyaEOIRwA682cp5PGcttktYhULRcYQ8Xk0I8X1gDnCNbKUeWs1aVCiim1AnKl0PPARcJaVsCI9ICoWiMwg1ZvACkAxsEELsFUL8MQwyKRSKTiDU8WqXtb2XQqHoCoRzNUGhUHRhlDJQKBSAUgYKhcKFUgYKhQJQykChULhQykChUABKGSgUChdKGSgUCkApA4VC4UIpA4VCAShloFAoXChloFAoAKUMFAqFC6UMFAoFoJSBQqFwoZSBQqEAlDJQKBQulDJQKBRA6OPVnnCNVtsrhPhECNE3XIIpFIrIEqpl8BspZZ6UMh9YDzwaukgKhaIzCHW8Wo3H20RAzUNQKLooIY9kF0I8CXwPqAamhyyRQqHoFEQrQ5C0HQIYr+ba72HALKV8zM959PFqwHDgaADy9QLKA9ivM4l2GaNdPoh+GaNdPghcxoFSygxfH7SpDAJFCDEQ+FBKOSYsJ9TOuVtKOSFc5+sIol3GaJcPol/GaJcPwiNjqKsJQz3ezgWOhHI+hULReYQaM3hKCDEccAKngCWhi6RQKDqDUMer3RIuQfywqoPPHw6iXcZolw+iX8Zolw/CIGPYYgYKhaJro9KRFQoFECXKQAhxvRDiqBDimBDiZz4+F0KI37k+3y+EKIgy+Ra75NovhNguhBgbSfkCkdFjv8uFEA4hxPxok08IcbUrtf2gEOLfkZQvEBmFEKlCiL8LIfa5ZLwjwvK9IoQ4J4Q44Ofz0O4TKWWnvgAjcBwYAsQB+4BRXvvMBv4BCGAS8FmUyfctoIfr529HUr5AZfTY71/AR8D8aJIPSAMOAQNc7zOj7W8I/Dfwa9fPGcAFIC6CMk4DCoADfj4P6T6JBstgInBMSnlCSmkD1gA3ee1zE/C61NgJpAkhsqJFPinldillpevtTiA7QrIFLKOLHwPvAuciKRyByXcbsE5KWQQgpYxGGSWQLIQQQBKaMrBHSkAp5WbXNf0R0n0SDcqgH3Da432xa1uw+3QUwV77LjTtHEnalFEI0Q+4GfhjBOVyE8jfcBjQQwjxqRBijxDiexGTTiMQGV8ARgJnga+A+6WUzsiIFxAh3Sch1yaEAeFjm/cSRyD7dBQBX1sIMR1NGUztUIl8XNrHNm8Zfws8JKV0aA+2iBKIfDHAeOAaIB7YIYTYKaUs7GjhXAQi43XAXmAGkANsEEJskc0L9jqTkO6TaFAGxUB/j/fZaJo32H06ioCuLYTIA14Gvi2lrIiQbG4CkXECsMalCHoBs4UQdinl+1EiXzFQLqWsB+qFEJuBsUCklEEgMt4BPCU1B/2YEOIbYATweWREbJPQ7pNIBmn8BD1igBPAYC4GbkZ77XMDzQMjn0eZfAOAY8C3ovVv6LX/a0Q2gBjI33Ak8E/XvgnAAWBMlMn4ErDc9XNv4AzQK8L/14PwH0AM6T7pdMtASmkXQiwF/g8tovuKlPKgEGKJ6/M/okW/Z6PdcA1oGjqa5HsUSAdedD157TKChS0BythpBCKflPKwEOJjYD9aevvLUkqfS2idJSPwBPCaEOIrtBvuISllxKoZhRCrgauBXkKIYuAxINZDvpDuE5WBqFAogOhYTVAoFFGAUgYKhQJQykChULhQykChUABKGSgUChdKGSgUCkApA4VC4UIpA4VCAcD/A1WVLN0nQDvBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# initialize likelihood and model\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = FirstGPModel(train_x, train_y, likelihood)\n",
    "\n",
    "# set to training mode and train\n",
    "model.train()\n",
    "likelihood.train()\n",
    "train(model, likelihood)\n",
    "\n",
    "# Get into evaluation (predictive posterior) mode and predict\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "observed_pred = predict(model, likelihood)\n",
    "# plot results\n",
    "plot(observed_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the kernel doesn't perform well. This is due to the lack of a lengthscale parameter, which we will add next."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding hyperparameters\n",
    "\n",
    "Althogh the `FirstSincKernel` can be used for defining a model, it lacks a parameter that controls the correlation length. This lengthscale will be implemented as a hyperparameter. See also the [tutorial on hyperparamaters](https://docs.gpytorch.ai/en/latest/examples/00_Basic_Usage/Hyperparameters.html), for information on raw vs. actual parameters.\n",
    "\n",
    "The parameter has to be registered, using the method `register_parameter()`, which `Kernel` inherits from `Module`. Similarly, we register constraints and priors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import positivity constraint\n",
    "from gpytorch.constraints import Positive\n",
    "\n",
    "class SincKernel(gpytorch.kernels.Kernel):\n",
    "    # the sinc kernel is stationary\n",
    "    is_stationary = True\n",
    "    \n",
    "    # We will register the parameter when initializing the kernel\n",
    "    def __init__(self, length_prior=None, length_constraint=None, **kwargs):\n",
    "        super().__init__(**kwargs)\n",
    "        \n",
    "        # register the raw parameter\n",
    "        self.register_parameter(\n",
    "            name='raw_length', parameter=torch.nn.Parameter(torch.zeros(*self.batch_shape, 1, 1))\n",
    "        )\n",
    "        \n",
    "        # set the parameter constraint to be positive, when nothing is specified\n",
    "        if length_constraint is None:\n",
    "            length_constraint = Positive()\n",
    "\n",
    "        # register the constraint\n",
    "        self.register_constraint(\"raw_length\", length_constraint)\n",
    "        \n",
    "        # set the parameter prior, see\n",
    "        # https://docs.gpytorch.ai/en/latest/module.html#gpytorch.Module.register_prior\n",
    "        if length_prior is not None:\n",
    "            self.register_prior(\n",
    "                \"length_prior\",\n",
    "                length_prior,\n",
    "                lambda m: m.length,\n",
    "                lambda m, v : m._set_length(v),\n",
    "            )\n",
    "\n",
    "    # now set up the 'actual' paramter\n",
    "    @property\n",
    "    def length(self):\n",
    "        # when accessing the parameter, apply the constraint transform\n",
    "        return self.raw_length_constraint.transform(self.raw_length)\n",
    "\n",
    "    @length.setter\n",
    "    def length(self, value):\n",
    "        return self._set_length(value)\n",
    "\n",
    "    def _set_length(self, value):\n",
    "        if not torch.is_tensor(value):\n",
    "            value = torch.as_tensor(value).to(self.raw_length)\n",
    "        # when setting the paramater, transform the actual value to a raw one by applying the inverse transform\n",
    "        self.initialize(raw_length=self.raw_length_constraint.inverse_transform(value))\n",
    "\n",
    "    # this is the kernel function\n",
    "    def forward(self, x1, x2, **params):\n",
    "        # apply lengthscale\n",
    "        x1_ = x1.div(self.length)\n",
    "        x2_ = x2.div(self.length)\n",
    "        # calculate the distance between inputs\n",
    "        diff = self.covar_dist(x1_, x2_, **params)\n",
    "        # prevent divide by 0 errors\n",
    "        diff.where(diff == 0, torch.as_tensor(1e-20))\n",
    "        # return sinc(diff) = sin(diff) / diff\n",
    "        return torch.sin(diff).div(diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define a new GPModel, train it and make predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAADGCAYAAADWg+V4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1t0lEQVR4nO2dd3hUVfr4P2cmk0x67yGUIC0kQOhFmgIKARusFF3FtQDLiqu4LmsBLGtdVlxWlJ99LWBB9AvqAorSl06o0kIJCekhPZnMnN8fk4wpk2QmmYQJnM/z5CFz59xz35lw33ve97xFSClRKBQKzZUWQKFQOAdKGSgUCkApA4VCUYlSBgqFAlDKQKFQVKKUgUKhABygDIQQeiHELiHEQSHEESHEYkcIplAoWhfR3DgDIYQAPKWUhUIIHbAVmCel3OkIARUKRevg0twJpFmbFFa+1FX+qEgmhaKN4RCfgRBCK4Q4AGQAG6SU/3PEvAqFovVo9soAQEppBHoLIfyAr4UQPaWUh6uPEUI8CDwI4Onp2bdbt26OuLRCobCDvXv3Zkkpg62912yfQZ0JhVgIFEkpX6tvTL9+/eSePXscel2FQtE4Qoi9Usp+1t5zxG5CcOWKACGEO3AjcLy58yoUitbFEWZCOPChEEKLWbl8LqVc64B5FQpFK+KI3YQkoI8DZFEoFFcQhzgQFVcvBoOBlJQUSktLr7QoCjvQ6/VERUWh0+lsPkcpA0WDpKSk4O3tTYcOHTDHlymcHSkl2dnZpKSk0LFjR5vPU7kJigYpLS0lMDBQKYI2hBCCwMBAu1dzShkoGkUpgrZHU/5mShkonJ6UlBRuueUWrrvuOmJiYpg3bx7l5eUAfPDBB8ydO/cKS1gXLy8vq8e1Wi29e/cmNjaWXr16sWTJEkwmU4NznT17lk8//bQlxKyBUgYKh5OWlsaIESO4dOlSs+eSUnL77bdz6623cvLkSU6cOEFhYSFPPvmkAyS1TkVFRYvN7e7uzoEDBzhy5AgbNmzgu+++Y/HihhN9W0sZIKVs9Z++fftKRdvg6NGjdp8ze/ZsqdFo5OzZs5t9/Y0bN8rrr7++xrHLly/LgIAAWVRUJN9//305adIkOW7cONmlSxe5aNEiKaWUhYWFcvz48TI+Pl7GxsbKlStXSiml3LNnjxw+fLhMSEiQY8eOlampqVJKKUeMGCEXLFgghw8fLhctWiTbt28vjUajlFLKoqIiGRUVJcvLy+WpU6fkuHHjZEJCghw2bJg8duyYlFLKM2fOyEGDBsl+/frJp556Snp6elr9PLWPnz59WgYEBEiTySSTk5PlsGHDZJ8+fWSfPn3ktm3bpJRSDhw4UPr4+MhevXrJJUuW1DuuNtb+dsAeWc99qZSBokHsUQZ6vV5izlit8aPX65t8/aVLl8pHHnmkzvHevXvLgwcPyvfff1+GhYXJrKwsWVxcLGNjY+Xu3bvll19+Ke+//37L+Ly8PFleXi4HDx4sMzIypJRSrly5Us6cOVNKaVYG1ZXXpEmT5E8//WQZ94c//EFKKeXo0aPliRMnpJRS7ty5U44aNUpKKeXEiRPlhx9+KKWUctmyZTYrAyml9PPzk5cuXZJFRUWypKRESinliRMnZNV9smnTJjlhwgTL+PrG1cZeZaDMBIXDOHPmDNOnT8fDwwMADw8PZsyYQXJycpPnlFJadYZVPz5mzBgCAwNxd3fn9ttvZ+vWrcTFxbFx40aeeOIJtmzZgq+vL7/++iuHDx9mzJgx9O7dm+eff56UlBTLnHfeeWeN31etWgXAypUrufPOOyksLGT79u1MmTKF3r1789BDD5GWlgbAtm3bmDZtGgB333233Z8RzDEdDzzwAHFxcUyZMoWjR49aHW/rOHtRcQYKhxEeHo6Pjw+lpaXo9XpKS0vx8fEhLCysyXPGxsby1Vdf1TiWn5/PhQsXiImJYe/evXWUhRCCLl26sHfvXr777jsWLFjA2LFjue2224iNjWXHjh1Wr+Xp6Wn5fdKkSSxYsICcnBz27t3L6NGjKSoqws/PjwMHDlg9vyke/DNnzqDVagkJCWHx4sWEhoZy8OBBTCYTer3e6jn//Oc/bRpnL2ploHAo6enpzJo1i507dzJr1qxmOxFvuOEGiouL+eijjwAwGo089thj3HvvvZYVyIYNG8jJyaGkpIQ1a9YwdOhQUlNT8fDw4K677mL+/Pns27ePrl27kpmZaVEGBoOBI0eOWL2ul5cXAwYMYN68eSQmJqLVavHx8aFjx4588cUXgPmJfvDgQQCGDh3KypUrAfjkk09s+myZmZnMmjWLuXPnIoTg8uXLhIeHo9Fo+M9//oPRaATA29ubgoICy3n1jWs29dkPLfmjfAZth6Y4EB3N+fPnZWJiouzcubPs1KmTnDt3riwtLZVSSvn+++/LKVOmyPHjx9dwIP7www8yLi5O9urVS/br10/u3r1bSinl/v375fXXXy/j4+Nljx495IoVK6SUZp9B1ZgqvvjiCwnIn3/+2XLszJkzcty4cTI+Pl52795dLl682HK8yoH44osv1usz0Gg0slevXrJHjx4yPj5evvrqqxZH5YkTJ2RcXJwcOHCg/Otf/2qZo7y8XI4ePVrGx8fLJUuW1DuuNvb6DBxez8AWVD2DtsOxY8fo3r37lRZD0QSs/e1atJ6BQqG4OlDKQKFQAEoZKBSKSpQyUCgUgFIGCoWiEkcURG0nhNgkhDhW2V5tniMEUygUrYsjVgYVwGNSyu7AIOCPQogeDphXoQDMkX3VQ3wrKioIDg4mMTHxCkp19dFsZSClTJNS7qv8vQA4BkQ2d16FogpPT08OHz5MSUkJYI44jIxU/8UcjUN9BkKIDpgrJav2agqHcvPNN7Nu3ToAPvvsM0tSEEBRURH33Xcf/fv3p0+fPnzzzTeAuQ7A9ddfT0JCAgkJCWzfvh2An3/+mZEjRzJ58mS6devGjBkzuBLBd86GwxKVhBBewFfAI1LKfCvvW9qrRUdHO+qyilbkkUegnhydJtO7N7z+euPjpk6dyrPPPktiYiJJSUncd999bNmyBYAXXniB0aNH895775GXl8eAAQO48cYbCQkJYcOGDej1ek6ePMm0adOoinzdv38/R44cISIigqFDh7Jt2zaGDRvm2A/XxnCIMqhsxf4V8ImUcrW1MVLKFcAKMIcjO+K6imuH+Ph4zp49y2effcb48eNrvLd+/Xq+/fZbXnvN3NGvtLSU8+fPExERwdy5czlw4ABarZYTJ05YzhkwYABRUVEA9O7dm7Nnzypl0NwJhDlv813gmJRySfNFUjgrtjzBW5JJkyYxf/58fv75Z7Kzsy3HpZR89dVXdO3atcb4RYsW1Zvq6+bmZvldq9W2aKmztoIjfAZDgbuB0UKIA5U/4xs7SaGwl/vuu49nnnmGuLi4GsfHjRvHv/71L4vdv3//fqAFU32vUhyxm7BVSimklPFSyt6VP985QjiFojpRUVHMm1c3jOXpp5/GYDAQHx9Pz549efrppwGYM2cOH374IYMGDeLEiRM1ipco6qJSmBUNolKY2y4qhVmhUDQJpQwUCgWglIFCoahEKQOFQgEoZaBQKCpRykChUABKGSjaCJcuXWLq1KnExMTQo0cPxo8fXyO82Fa2bNlCbGwsvXv35uLFi0yePNnquJEjR3KtbX+rjkoKu/jnBvtvwIb485gujY6RUnLbbbdxzz33WBqVHDhwgPT0dLp0afz86nzyySfMnz+fmTNnAvDll1/aL/RViloZKJyeTZs2odPpmDVrluVY7969GTZsGI8//jg9e/YkLi7O0huxvhTld955h88//5xnn32WGTNmcPbsWXr27AlASUkJU6dOJT4+njvvvNNSOwHMiVCDBw8mISGBKVOmUFhYCECHDh1YuHAhCQkJxMXFcfz4cQAKCwuZOXMmcXFxxMfHW9rD1TePs6BWBi1AhdFETnE5uUUGCssqKDMYKa0wUmowISVoNQIXjUCrFXi6uuDj7oKvuw4/d1fcXbVXWnyn4/Dhw/Tt27fO8dWrV3PgwAEOHjxIVlYW/fv3Z/jw4YD1FOX777+frVu3kpiYyOTJkzl79qxlruXLl+Ph4UFSUhJJSUkkJCQAkJWVxfPPP8/GjRvx9PTk5ZdfZsmSJTzzzDMABAUFsW/fPt58801ee+013nnnHZ577jl8fX05dOgQALm5uY3O4wwoZeAAcorKScktJiW3hIz8Ui6XVGBqYpi3t96FcF93Ivz0RPq7E+zl1qSGntcCW7duZdq0aWi1WkJDQxkxYgS7d+/Gx8fH7hTlzZs38/DDDwPmdOn4+HgAdu7cydGjRxk6dCgA5eXlDB482HLe7bffDkDfvn1Zvdqcvb9x40aLOQPg7+/P2rVrG5zHGVDKoAmYTJJzOcWcTC/gXHYxhWWOS38tKK2goLSAE+nmRptebi50DvWiS6g3Eb56i2JIS0tj6tSprFq1qlldjtsCsbGxVm37hvJqmpKiXF/r9zFjxvDZZ581eJ3q15BW2sg3No8zoHwGdnAhp5iNR9NZseUMa/Zf5EhqvkMVgTUKyyo4cD6Pz3df4N2tyWw/lUVhWQXPPfccW7du5dlnn23R6zsDo0ePpqysjP/3//6f5dju3bvx9/dn1apVGI1GMjMz2bx5MwMGDGjSNYYPH27pnnz48GGSkpIAGDRoENu2bePUqVMAFBcXN7qLMXbsWJYtW2Z5nZub26R5WhulDBqhrMLIgQt5fLTjLO+t38dDUyeSXqvNeH52Bsseu4v8nMwWlaWgtIIRsVF463UsX74ck8nE8uXLEULg7u7eote+kggh+Prrr9mwYQMxMTHExsayaNEipk+fTnx8PL169WL06NG88sorTV4lzZ49m8LCQuLj43nllVcsSiU4OJgPPviAadOmER8fz6BBgyyOwvp46qmnyM3NpWfPnvTq1YtNmzY1aZ7WRqUw18PlEgP7zuVyNC2f8goTAF++sYgd61YyeMJUJj+8yDK29vH87Aw++vuj/P7Jf+ITEOxQufKzM/h2xcsc2r4RQ1kprm56Eifdwr/feL1FzAWVwtx2sTeFWfkMapFTVM6u5Bx+vVRgcQL+JTGeivIyy5jtaz9j+9q6tl/VcaHRgJSs//jfNZSGI/AJDEHv4UVFeRkurm4Yysu4VKJhV7qJIV5lBHm5NT6JQmEFZSZUkltUzrqkND7acZZjafk1dgOe+nAjCaMS0bmZa+jp3PQkjJ7IY8vX1DhehTSZkFKyfe1nPDq2K39JjHeorAV52QxJnMa8pZ8zJHEaBblZnM4o5OOd59h4NJ1SgyrvpbCfa35lUFhWwc7T2RxJza+zHVh9uV/9aVxRXobew4vImO51jgNoNFpMJiM6Nz1xQ8cw6cEnGpXDHtNi5sLfnFN3/Gmh5Xcp4dDFy5zKLGRY5yBiI3zUtqTCZhyyMhBCvCeEyBBCHHbEfK1BWYWRrSez+GBbMocuXrYaF7D+kzdJPryH9R//2+rTGMxPaYSoYUaYTOYns6GsFL2HV52b25rDsepaa995jWWP3cXF08ea7JQsKTey4Wg6X+xJIbuwrPETGkE1GGl7NOVv5hAHohBiOFAIfCSl7NnY+CvpQJRSciQ1n+2nsygqs76cru0jqMLF1Y1X1ibVOV7bqSc0Wrr2HYqXXyBlxYU1nuRQ0+G4a/1qq9cCsxe9trPSXlw0giGdgwjXlTBt2jS74xKSk5Px9vYmMDBQrTLaCFJKsrOzKSgooGPHjjXea3EHopRyc2VrNafmYl4JP/+aQUZ+w0/Lpz7cWOPmdnF1Q+/uyUMvv19jXH1mhNFQTkBoZJ2b2JojEkAIDVKa6shR5XfYvvazehVRbWqbGxUmyeYTmfzw9m9xCW+++Waj81QRFRVFSkoKmZktu22qcCx6vd4SgWkrDttarFQGa+tbGdRqr9b33LlzDrmuLRSVVbDlZCbH0gpsPufLpQvZ8d0qtDpXyw08JHFavVuKBblZePuH0K3/TLavPUhuuh+RnSdTkONCcYGG4gItRfmCivIKTKZioBQox8OnDDf3dHLTtwDngCRgH2AAzKuR+GFjGTXlD3z95gt1fAq1b/4qmfrecAs56Rc5f/wgFYbyOp9Pr9fXSMZRXBs0tDJoNWVQndYyE0wmyYGUPHaczrbECtjK+4vncnjHj0hTQ+cJIB4YXfkzHPCxvOsXbMAvuAIPb2Plj4lTB7aRmnwOoXFHmnT4BvWipDAQQ1koUlYt1Eow967dBmxk8IQIhJANxjgghHVZhSBh5ATLKkfnpuf6G8fzwdv/ol1khF3fiaLtc00qg9S8En48nkFWQdMdaLV9ATo3PT2HjOO63vP5aVUBWal9gRAA3DxSiBuipVOcIKxDOWHR5ZSXXqqzQ/D+4rn4BAQzaPyd7PxuFfk5mcxcuAyjEd556jmgL5fOh1F0OZaK8u6YLblLwJfAKswKoml/syoTZvCEqcx68kUS48Lx93Rt8vejaHtcU8qgpNzIlpOZHE3Lx9aPVrXUvm3Ok3WW4lXmgsYlHqPhXnRu92Ao80ejLcVk/AaNdj0m43qGJF5fx0dQX8SirZQWC/ZvqmDDp9nkZSYAHsA5Ijv/yK2zwtixbnkNReUbGEpW6rkaPgidmx4Pbz+u6zOYEbffW0MBubpoGNMjlC6h3nbLpmibtLgyEEJ8BowEgoB0YKGU8t36xreEMqjaJdh6KouScvuCbqpu2pDoGDLOn7bY2zOeeJ0Pn/uRvKw7yM/ugdAY8Pb7H7fM6sS+n2bhF+xd5wkP9u9GNCrf0oVsX7cWjfZWTMb7gFG4upnwD/2R9PNzcXE9h9FQjn9oJN37D6cgL4ekLT8ghAaQjSqj/h0CGNpZ7RZcC7TKysAeHK0MMvJL+el4BmmXS+06r76bFjyBh4BHgUj8Qw0MScxj4Lh8vPyMjQYIWTMvqoKPmpKrUN202Lz6Aw5tz6Rb/y85+Is/UroQOzgFd4/XKCs9zsyFy+o1RRqiY5AnN/UMQ69TxVWuZq7a9mqlBiM/Hkvn013n7VIEVUE/85auImFUIi6uVfH8vsBTmL36/wCOA4nk5/hyw525ePmZVxzVg5GsUTt/oCpiESlrBBLZmu04c+Ey7vjTQiJjuqFz01NWvA0Prz+zaOUFRtyRx697oti/+XX8QlZRkKutMf6OPy1sVBEAJGcVsXLXeYcEKSnaJm1yZWAySQ6nXmb76Wy7TQKoacsjJdvXrQXxZ5CPYVYI/we8gM7tYI0nuj3Lf2tPZ2//oBo+BHt8Cg1de8F7R1n/cSC71vvgpjdx873ZDJ2Yh6YJD3k3nYYJceG0D1Qdi69Grioz4Xx2Mb+czGzSLkHdG0qHOfThaSAUN/f1uLi+SNHln63a201d/v9lQpzVvf7aWFMq1Z2bmz5/p8Frp5/XsWZ5CL/u9SSqcymT56UT3dX+70kjBCO7BtOrnZ/d5yqcm6vCTMguLOPbg6l8tS+lyduFVdmHLq564A7MZsAyOvTQM2/peV78pgOdeuoYOnE6j775dY0cBKh/+d+YH6D3iPGAOYEJzDe9f2ikxTypyoJ86qMf65xbZZLsWLeq0WuHRht48O8X+f2TqeTnuLD04Wi++lcIZSX2OQZNUvLT8Qw2/ZqByaTyEq4VnD5r8XKJgZ1nsjmeVtBgkVFbsv58AkMwGmKpKJ8PjACS6D7gJe5/7naqHOn1ZQRWUZWwVH35Xx+1VyJVCUwV5WW4urljNJTXe2PXG7qs0TBv6ef1XlsI6D2ikG79ivn+o0C2rvHj+B4Ppv/lEh1j7XOwHjifR36JgZt7huPq0maeG4om4tRmwuYTmRy4kIfRhqdTQ/Z3fnYG7z/7LL6BK0jaGo2LrpCRU5Ipzn+VgtwMmxxsTaGhBKYTe7cSN3RMvR5/R+1InDnkzqevhpKbrmPU73K56e5sXFzt+5uH+LhxS+9IvNyc/tlxTbL1ZJbNW8NtttLR8Uv5jSqC+qoQVdnfJhN89PfjnDu2EvAG3iBh1BHG3/s4+dmz+Ojvj5Kfk+nw8mRQ16yoL4HJ2grEXpOkvpVRp7gS5r91jm/fDuanVQEc2+3J7/+WSmi0webPkZFfxspd57m1T6SqpOREyEpzLinlMkM7BzZ7vja/9quvCtFTH/3I4+PvYf5NuZw5dD/mBKBewJ/Ztf4dHh3blcUzRjS4RdhU3F21RPq706udLzpDAVPuvo+vvv+Ju2bej48s4rY+kUzqHcGE+HBGdQshob0/MSFeBHm74aL5TbvXV0PBGg1td+o9JL/7cwb3P3eR/Gwt/5zbnr0/2Rd1WFBawed7LnAhp9iu8xQtg9Ek+f7wJZJSLjtsTqc2E1ZsPl1vzYHqVM8wNBrK6T92JvAcu/4bhosuD5PpUUzGDxudpykRgkJAkJcb0QEetAvwIMTbDc9mLKcrjCbSC8q4mFtCal4JF/NKGkyyamy7s/aKIS/Lhf/8PYzkwx4MGp/HbXMy0dlhNmg1gjE9Quke7tP4YEWLYDCaWJeURnJWkeXYIzded3WbCdawthyueoIOvPlOvnnrGLvX342UYcBy/EPfJeviAYRGgzSZLP8GRbTncnZ6HXvcFlw0gvZBnnQJ9aJ9gKdDW6K5aDVE+rkT6WcufW4wmjibVcSJ9EKSswoxGGveuLVrL9T+LNVXDJMfXoRfUAVzXk3h+w+C+GlVAGePuKDT38UfFj9uk6lkNEl+OHyJgtIKBnQMcNjnVthGqcHINwcukppnnzPYFtqcmWBtOTxz4TJG3/kc//1oNKeTnkDKDGAQMJfMlP1IKZFSEtb+OmLiBzB04nRMpgq7tgiFgOgAD8bGhvLA8E5M6hVBtzCfFu+NqNNquC7Umwnx4Tw0IoaxsaEEe/9mt9fnW3j+9zfw6NiubF/7WZ3irFotJP4hi/ufu0hmquDCr+/wxeubbJYpPzuDabfcxNurf2L48BFcquwjkZaWxogRv71WOJaCUgOf77nQIooA2tDKoD5HoVbnzsT7z7JmuQfm+gKPAUuB38yLhNET63jh3188l+79RzS6Reim09Aj3IdeUX5XPN1Xp9UQG+FLbIQvF3KKOXAhj9OZhVa3OxtbMfwWCBUDfMORnX/l0bHz0bgsp0P3Xg1u0VYp5IWPPkTG+dMsXLSYt99aXqPLkz3VlBSNk1NUzup9KRSUtlwHrzbjM7C21RYT/yD52c+Tesabzr1zcNM/zon9n2IoM2vOqv4FTUkh9nXX0a+DP93CfJx6jz27sIztp7M5lVG3vXdtX0r17+HTV55gz8Y1lZWc3RGaT5CmSQSGbyI7bTxDEm9rtGxbY6hqSo4hNa+Ebw+mNhh67wifgfP+L6+kKpkHISzLYa3OH0PZsxzfs4T8HHfuXpDG7Jez8A0spqK8rDKUGOKHjW3UC1+bQC9XbuoZxr1DOhAf5efUigAg0MuNib0imDYgmugAjxrvWduN+MuEOB4d25U9G9cAVYFQhUjTrcAzZKeNAjayfe1GHh3blccnxFnm+y2Cs+b2olbnSmBYJHq9eUfHw8ODGTNmkJyc3HIf/BrhVEYBq/elNCkHx16c3kyoXa68+4AXSTk1l/xsT7QuHzHnta6ERfsD1qMDre3hWyPA05UhMYF0DvFqk3n9Yb567ugbRXJWEZuOZ3C5xGA1mrLmisBoCYTKSEkmJ+05hDiOlB9gLruWSJ+RXSxzVPknjIbyGg5ZU4UBF1d3ysrKcHPTU1paio+Pz1XfHbql2Xc+l80nMhst0pOfncHIkQ80uyO305oJ7u7ulJZWd5SEAq8DUwlpV0Zw1Msc3bmo2aXEvfUuDOoUSI9wHzSatqcErGEwmtidnMOec7mWoC37lvj9gG8xV1b6HbDeslVZlY2ZfuEMhbnZePkHEtquE4e2bSBu6BiGJt7Jue3fUp6fw+rVq1vmA17lSCn55UQm+8/n1TheX2DZl28sYud3q3jooYca9dW0yazFAwcOcP2oGygpKsFomAm8ALiD+DvIl4CaWYBanSuvrjtkswyuLhoGdAygTzs/XLTObQo0lZyicjYeS+dibkm9odGu7p6cO7qf4oI8s69FaPD2D6Qg1wPkGiCOdl3f5g+Le9kVpTmgYwBDYlT1JHsprzDxw5FLnLbmA3qj+ZWv26TPYMWKFRTmdcRo+AX4N7AHiAP5LNUVgajMBOwzcrxN8woBsRE+3DukA/07BFy1igDMps+UvlEM7xJEQHBojS1IpImA0Ejufep1YgeOtBwXSDy8fBGcR6u7AfieC7/O4ZfV3cnLtL31/K7kHL47dAmD0b6q1NcyVVuHtRXBXxLja2wT79m4hjOHdlNRYSBhVCJa3XhgYrN9NY5qr3aTEOJXIcQpIcRfmzOXu7s7QgiWL48GdgFRwFRgDHASnZueoIj2lvGyMhNwz4Y1jTY5jfDTM21ANGNjw5oVJdiWEELQt30A0wZGU16YW2+LuOrHSwovMyRxGo+88T79bvgAF927bPo8gH8/ruXMoUM2h2+fSC/gy70pFJa13HbY1cKly6Ws3HWBTCvp+VWO2zpIyb5NURgN3yDEk5SUNM9X02wzQQihBU5gvltTgN3ANCnl0frOachMSEtLY/78+XzxRQQGQzQurs/jF+ROdtp5yxaZf2gkMXH9KcjL5sS+bZiMxgaz+txdtaoRKebowR2ns9lzLsfmytFfvrGoMn16PvAqsAW4FcixOXzby82FxF7hhPu6N1X0q5rjl/LZeDS9TnRpdcxFcVdWq3ytRaNdisn4R3yD9rBujSuffvo2aWlpDfpqWjoceQBwSkp5pvJiK4FbgHqVQUOEh4fj4+OD0bgEnauOCkM5JlNknV2CaY+/xJdLFyJNpnqjCM0mgS/XXxekCn1izisYdl0Qkf7u/PfIpQa3q+o6HF8DzgL/AbGd2EGvMWXevTZdt7Csgi/3pDC6ewixEb7N+ARXFyaTZPPJuo5CaxTkZTN04vTKytfbgFWYjDfjql/OvDfiGTKkL0OHNi/hzhErg8nATVLK+ytf3w0MlFLOrTXO5vZqt99+O+Hh4QT0vZlNaz6rt7pvQ1WAA71cubF7KBF+6mlkjfxSA98fSqs3tNWaw1GajGi1IzEav0LnqmXukjzadbGv6lSfaD+GXxd81ezcNJXi8grWJaWRkmtfUNbbC57lwokXKS2KxMvvGfJzXmBI4jS2fvtJs4OOHKEMpgDjaimDAVLKP9V3jqOzFquj1Qj6dwhgQMcAtNf4f7jGMJkkW05lse9crtX3a/ebDG3fmbv++g9+WrWZwzvmIwjh90+l0WNgkdXz6yPCT8/4uHC89TpHfIw2R9rlEtYlpdkdWnzumJ53F0ZQmFeIuWzfTzXetyXis6V3E1KAdtVeRwGpDpjXbiL89MwYGM3gmEClCGxAoxGM6BLMTT3D0Gnrfl/VHYtDJ04nJKojkTHduPtvD/Lkh0WERJfz7sIIdqyzb+mfmlfKJ/87z67Dp6+pxCYpJbvP5vD57oZzDGqX0M/PzuCl+9/j3/MjcdWb+ONrp0kY5VGjhsf06dObHfHpCGWwG7hOCNFRCOGK2fX/rQPmtRmdVjCiazC/69eOQFWJx266h/vwu37t8NbXdCE11H/BJ8DIH1+7QPf+RXyxNJS17wbRYI/aWpSUG3n0b0+zZctWFi9e7KiP4nRUZXKeOZfC1/svsvVklqWWZ319M6pH3UoJ7y2+QMb5l3D3OsUjb1wgJt67TqaqIyI+HdVebTzm8EAt8J6U8oWGxjvSTIj0d2dsj1D8PFQD0eZSXF7B2qQ0LtphxxqN8PW/Q9i+1o8+I/OZNj+90RqL9UVDXo2JTXPmzOGtt9/m+olTufWPNUPjGw8i0gFvAzOBzyr/NSuA7v2H1/CXhbiU2BTx2SYjEKFhZaDTCoZ2DqJ3O79rervQ0RhNko3H0jmamm/zOVLCpi/8WftOMJ3iipm5MBVPn/qXCdYyUHsNG8NLL7/KDX27OuJjXHHqhtObqUryshoaLgQJIydUfi+eIL4GeT1C8zzS9HSD2+fXRNaiNcy+gfb0ifZXisDBaDWCcbFhDLsuCFu/WiFg9O9yuXtBGueO63njkWgyL/7mHKy9HLZWkMXV3YukHMGa/RfJL7W9WKuzsmn3IQbcOMlqbc6Gg4jWYiiLBv4Hsj8+gY+CfMauPh1NpU0pg6p98il9213xQiNXO/07BJAYH27VsVgffUYVMPvlFIoLNCx9OJrTSeZt3SobeO07r7Hssbu4ePoYB7b8QN8bb60TDZmcVcR/dpxjz9mcNtnApbCsgh8Op7EjzYSLm4fValpVyhCwpNuDWWF4+NyJ1mUv7l5RxA19CeTHNhfFbS5txkwI8nbjptiwGiW/FC1Pen4p3xy42KDvpnY23dmj+Sx/IghDWQRwP/CfOucIIRrNOA3ycmV091BLPUhnpsJoYv+FPHYl51gK2DYUB1P1njmI6AfMz+WHgX8Q3tHAfYsuonNNbbQxUBWOMBOcXhkUlxvpE+3P0JjAqzqpyJm5XGLgmwMXyS603i+ydgMbcwjzD3j7/0xBbjzwIubu1tb9CA2FNQsBnUO8GNwp0Cl3iqSUnM4sZMvJLPKK7Tdv3l88Fy/fcArynuXw9vb4BGxjwfshuLlLuxrzXvXK4NP/nef664JoV6uCj6L1KTUY+b+DqTUi5hqvkeACLAMeAtYB04HfHJMurm7EDxtrU5cojRB0C/dmUKdAfN2vfLCSlJJTGYX8LznHanKRreRna/nw+QiSj7gzZno2436fzV8n2d7tG8y5Hw8M72TT9dqsA3FKvyilCJwEvU7L7QlRdA//rflK7QY2dRvKuuAbtIjobm+BGIe5etJ1gLk+pdFQbpNDLD87gzcencHOw6f5cPtZfjh8ibTLV2YLssJo4lhaPh/vPMfapLRmKYLkI3r+8cf2XDzlxt1/S+Xme7PRaBpuDFSb6AAPZgyKbrIM1XHqPF6dMguciqqdBm+9jl3JOVbbx9VuKBs7aBSTHx7Nskef4PyJhZgq9uIT+BhBkYcIbdfJptoItXs/HEvL51haPqE+enq186VLqHeL/1/JLCjjcOpljqcVUGpoXj1CKWHrN35883YwAaEGHnoxhYiOv5lgtrTWEwIGdQpkYMcAh+2oObUyUDgfQpjjO7z1Lmw6nlmn7uShbRusdqmeu+QhctKzeH9xBBdPraDvDdncdE82WivJpFUOydpBOLX7aKbnl7L+SCmbjmfQLsCDmGAvYoK9HNLLwmSSXMov5Vx2MclZRaTnO6ZXQVmJ4Iuloez7yYfYQYVM/8sl3L3q+lIa6vbt4aplfFy4w1fNTu0zUDg3ZzIL+f7wpQbbv9XGUC74+s1gdn7nR0x8MXcvSMMnsOaTtnpknslYYVcnao0QBHm7EuKtJ8TbjRAfN/w9XBtMYTeZJHklBnKKysguLCezsIwLOSXNXgHUJuWkG//5ezhZaTpuujubG6bloKm2oKmvxmF1wnz1JMY3PcnrqmqvpnAeOgV7MblvVKNbj9XRuUp+90gGnXqW8OXSUF6b3Z4ZT1yia9/iOg7JqnLugM1BNyYpycgvIyO/pi2v1Qj0Og16nRadVoPBaMJglBiMJsorTI12+7YVaze0lLDhExd++E843v6S2a+k0Dm+rs+jtjlUm9gIH0Z3C2mxXTVllCuaRaiPnjv7RxNgZxBYvxsLeGTZeTx9jby9IIrV/w5mzitf4OkbUM0Bqcc3KIx+Y25rdtCN0SQpKjOSXVjOpculZBeWk19ioKTc6DBFAHXb/+XnaHn3mQh++KgTyB/o1u9PhESeqxGRWbvGYfVWeGBWZKO7hTA2NqxFt9eVmaBwCKUGI98eTLUryQmgvEyw7t0gtqzxx93rIiWFtwB7LQ7J5pbCby2sb7POAN7AXHL+cczbrL8xJHEakx9eZDVXo8ocCgsLY0J8OFH+jvEPtNmtRUXbQa/TckdCFN3CvBsfXA1XN8mO78KBGykplMAO4FkqygUI0aLht46gKu9i3tJV1bYDwxFiLfAxoe0N9Bj4CDq3d+qcW7UCeP6eG63uHnTuEMW0gdEOUwSNoZSBwmFoNYKbeoYxsJN9rdrN++ruuLj2B1YCT+Pmfobpf0myWu7OmagyC3asW4Wr3hdD2YPAUaQcBfyZTj0fxS/ocmVbQLMppaks7189fqB2hWpZksfv+rXDpxWrQSlloHAoQgiGxAQxLjbM5mpTv7VtS8fF9QHgRoRG8slLsXz4XDi5GQ37ua0VCamvcIitNHZ+XTs/g53fL8DcAXw30At4nR3rPmX7upUgBI+88QVh7TtjMhnrOESrCslEde7GkqX/YvvG71o9zkYpA0WL0CPCh6ERGpbPt+2GrFliLZSY+AcYNeUcSVt1/H1mB755O4jCPOvbg9WddlU38dr3/lHDkVeFrUqitiOwNpYmtLo+wFfARlz1Adwy63/0GfkGOrcU4Len/8JPNxMZ043gqI4MnTjdqkPUTafhlt6RDOho38rKUSgHoqLFmDNnDm+//TYjb5lO4uyn7T7fnPC0g+CoD8hKHY6rm2T47bmMvCMXdy+TXf0jqwKVGkv+qW9OjYuODt17W7YMf91bxLvPpFBhmAgUAC8y6OaL/O7PT9YoJGurEzTA05WJvSLs3pWxlzabqKRomzRU5ceWpivWb8huCPE8Ut6Bm7uJfjfm02v4aXZ+t8jiha+PhNETSdryX6t9Catu8tvmPMnXb77AmBmz+fjF+ZSVFFFRXmYpER8c1YHMlAt0H/A8Gs1cjuz0Ai7j6fsx9zzViYObP7akKDeUumyNTsGe3NQzDDeXlu/t0WLKoLJM+iKgO+by6Dbd4UoZXN1UdcVas2YNxcXFeHh4cOuttzJ17t84Vdi4Q6y+rbZRU/7Ayn98QVDEPziyI5gKgwatyxaMFcvRuKzHVJFtaTUP5mQopGTwhKmMnTHHah+I4KgOZF08R0h0DBnnTxMSHUP6uVO1JIrGXH9wJtAeyMZc8vNfwGXA3Pi3fbdeNtUeqEIIGNwpkAEOzC9o/Jott7V4GLgd2NzMeRRXEVVdsUpLS9Hr9ZSWluLr68vEwbFMiA/H1aXmfztbyqLpPbzYsW4Vqac/xdPnzzzzSTLR3T7DWBENfAoyA1f9L5hMDwAdAYgfNtZil1fNWbWCqOrRmZlyFikl6edOWf41FyIdjrnz9x7gHPAM8CtCOwO/kP64uL4KXLb4BPqMHN+gj6E27q5abusTycBOztOp2lHVkX8G5quVgaKKqq5YDz74ICtWrKjRAzCnqJx1SalkVRZLsWbHV19qL/njbch667BrgEGY+z/eSlWKtIuuAL3nCQaM60hodBlevkZ++vwlPLyNlBSlcurgIZDugCcQBHQDelT+xAHeQAWwA8QPID9GaFJASkLaxZBx4bSluYw1GjKJwnz1TIgPb9Vtwypa3GdgizKwp72a4urHYDTh5eFBuQ1FPGqbDS6ubnj7B1GQm0VFeZnFjJj4wBOUFERy+pA7Kaf0XDzpRmqyKyajrQvgDOAY5gXvj/iFHCd2YB/SL5yhMDcbL/9AQtt14tC2DcQNHcOg8XeyefUHnNi/g+KCvEaTqYSAvu39GRITdMWa/DQrUUkIsRGw1p3hSSnlN7YKIaVcAawA88rA1vMUVyc6rYazZ5O5f87DbPh+bZ0bqTq21E3Qe3jhGxiMb2A5YR3KqbLlP3/9OXZ+t4u4oVMZMO5BivK1/LL6Szy9dRTkXaCsOIPS4gv4+OdQVHASb78gy02fnxPFHX9aWEf26seqGgDv+G5Vg8lUnm5axsWG0T7Q0/FfpoNoVBlIKW9sDUEU1x7h4eFEhwVhNJSjayQrsXp+/+bVH7B30//R98ZbGXH7vXXy/aHujsShbc9yaNuzNu9o2ENDtQcAOgZ5MjY2FA9X504Sdm7pFFc96enpzJo1i/vvf4AX/vEvTp67YHVc9a05nZseaTTi6qa3tH6rzVMfbqw3+cfRVJetuixuOg0jugS3mTb0zVIGQojbMO+vBAPrhBAHpJTjHCKZ4pqgekuwLz9+l/T8UjYcTbdaW7D207525aPq2FI6rCXpFOzJDd1D8XJrO8/bZm0tSim/llJGSSndpJShShEomkuoj57pA6IZ0TW4zhakPYVCgTrJP47IgGwsnNlb78LNcWHc0juyTSkCUGaCwgnRaAQJ0f50CfXml18zOZFeANj/tK9v+d4c6qtG5KIRJLT3p3+HgDpKrK2glIHCafFyc2FCfDh98vzYdiqLlNySRp11LUVDJsrXu89wfedgfD2ufD+H5qByExRthvPZxWw7ncWly/ZXKral2Ghj59d2SA69YTzL3lhCbEx7u+e7UqhKR4qrguhAD6YNiGZy3yg6h3ihsSOMt7GU5MaobqK4upljHbq3D21TiqAxlJmgaHO0C/CgXYAH+aUGDqVc5lhaPgWlFVbH2rMD0RBhvnr0xkIeePAh5syeZQmxvppQZoLiqiA9v5TTGYWcziy05DxA/RmQjfV3dNEIwv3c6RDoQYcgT4KcsOlrU1B9ExRXPaE+ekJ99AzpHERRWQUZBWVk5JeSHuLFL76+De5AuOk0BHi44u/pSqCnK0FebkT6u19z7f2UMlBcdXi6udDRzYWOQeY8gA9EMbNnz+aBBx5g+Vtvk5qWxkMjOuGi0aDTCqdJIb7SKDNBobiGULsJCoWiUZQyUCgUgFIGCoWiEqUMFAoFoJSBQqGoRCkDhUIBKGWgUCgqUcpAoVAAShkoFIpKmqUMhBCvCiGOCyGShBBfCyH8HCSXQqFoZZq7MtgA9JRSxgMngAXNF0mhUFwJmlsQdb2UsiqRfCcQ1XyRFArFlcCRPoP7gO8dOJ9CoWhFHNJeTQjxJOYulZ80ME/1XotNElahULQczW6vJoS4B0gEbpAN5EOrXosKhXPT3I5KNwFPACOklMWOEUmhUFwJmuszWIa5kf0GIcQBIcRbDpBJoVBcAZq1MpBSdnaUIAqF4sqiIhAVCgWglIFCoahEKQOFQgEoZaBQKCpRykChUABKGSgUikqUMlAoFIBSBgqFohKlDBQKBaCUgUKhqEQpA4VCAShloFAoKlHKQKFQAEoZKBSKSpQyUCgUgFIGCoWiEqUMFAoFoJSBQqGopLnt1Z6rbK12QAixXggR4SjBFApF69LclcGrUsp4KWVvYC3wTPNFUigUV4LmtlfLr/bSE1D9EBSKNkqzqiMDCCFeAH4PXAZGNVsihUJxRRANNEEyD7ChvVrluAWAXkq5sJ55LO3VgK7ArzbIFwRk2TDuSuLsMjq7fOD8Mjq7fGC7jO2llMHW3mhUGdiKEKI9sE5K2dMhE5rn3COl7Oeo+VoCZ5fR2eUD55fR2eUDx8jY3N2E66q9nAQcb858CoXiytFcn8FLQoiugAk4B8xqvkgKheJK0Nz2anc4SpB6WNHC8zsCZ5fR2eUD55fR2eUDB8joMJ+BQqFo26hwZIVCATiJMhBC3CSE+FUIcUoI8Vcr7wshxBuV7ycJIRKcTL4ZlXIlCSG2CyF6taZ8tshYbVx/IYRRCDHZ2eQTQoysDG0/IoT4pTXls0VGIYSvEOL/hBAHK2Wc2cryvSeEyBBCHK7n/ebdJ1LKK/oDaIHTQCfAFTgI9Kg1ZjzwPSCAQcD/nEy+IYB/5e83t6Z8tspYbdxPwHfAZGeSD/ADjgLRla9DnO07BP4GvFz5ezCQA7i2oozDgQTgcD3vN+s+cYaVwQDglJTyjJSyHFgJ3FJrzC3AR9LMTsBPCBHuLPJJKbdLKXMrX+4EolpJNptlrORPwFdARmsKh23yTQdWSynPA0gpnVFGCXgLIQTghVkZVLSWgFLKzZXXrI9m3SfOoAwigQvVXqdUHrN3TEth77X/gFk7tyaNyiiEiARuA95qRbmqsOU77AL4CyF+FkLsFUL8vtWkM2OLjMuA7kAqcAiYJ6U0tY54NtGs+6TZuQkOQFg5VnuLw5YxLYXN1xZCjMKsDIa1qERWLm3lWG0ZXweekFIazQ+2VsUW+VyAvsANgDuwQwixU0p5oqWFq8QWGccBB4DRQAywQQixRdZM2LuSNOs+cQZlkAK0q/Y6CrPmtXdMS2HTtYUQ8cA7wM1SyuxWkq0KW2TsB6ysVARBwHghRIWUco2TyJcCZEkpi4AiIcRmoBfQWsrAFhlnAi9Js4F+SgiRDHQDdrWOiI3SvPukNZ009Tg9XIAzQEd+c9zE1hozgZqOkV1OJl80cAoY4qzfYa3xH9C6DkRbvsPuwI+VYz2Aw0BPJ5NxObCo8vdQ4CIQ1Mp/6w7U70Bs1n1yxVcGUsoKIcRc4L+YPbrvSSmPCCFmVb7/Fmbv93jMN1wxZg3tTPI9AwQCb1Y+eStkKya22CjjFcMW+aSUx4QQPwBJmMPb35FSWt1Cu1IyAs8BHwghDmG+4Z6QUrZaNqMQ4jNgJBAkhEgBFgK6avI16z5REYgKhQJwjt0EhULhBChloFAoAKUMFApFJUoZKBQKQCkDhUJRiVIGCoUCUMpAoVBUopSBQqEA4P8DrGwjOCXQPogAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Use the simplest form of GP model, exact inference\n",
    "class SincGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super().__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = SincKernel()\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "# initialize the new model\n",
    "model = SincGPModel(train_x, train_y, likelihood)\n",
    "\n",
    "# set to training mode and train\n",
    "model.train()\n",
    "likelihood.train()\n",
    "train(model, likelihood)\n",
    "\n",
    "# Get into evaluation (predictive posterior) mode and predict\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "observed_pred = predict(model, likelihood)\n",
    "# plot results\n",
    "plot(observed_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because many kernels use a lengthscale, there is actually a simpler way to implement it, by using the `has_lengthscale` attribute from `Kernel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAADGCAYAAADWg+V4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1xklEQVR4nO2dd3hUVdrAf2cmk0x676FGakgIARGkCEhRDLgiFsCyutYVxcKu8lnAttZlFxddl1WxA4qKiq4LuCBdAQldIgQICSW9l2nn+2OSmDKTTDKTMIHze548ZO6ce+6bGe573/OetwgpJQqFQqE51wIoFAr3QCkDhUIBKGWgUChqUMpAoVAAShkoFIoalDJQKBSAC5SBEEIvhPhJCLFHCHFACPG0KwRTKBQdi3A2zkAIIQBfKWWZEEIHbAbmSCm3u0JAhULRMXg4O4G0apOympe6mh8VyaRQdDJc4jMQQmiFEGlADrBWSvmjK+ZVKBQdh9OWAYCU0gwkCyGCgC+EEAOklPvrjxFC3AXcBeDr6zu4b9++rri0QqFoBbt27cqTUobbes9pn0GTCYWYD5RLKV+1N2bIkCFy586dLr2uQqFoGSHELinlEFvvuWI3IbzGIkAI4Q2MB35xdl6FQtGxuGKZEA28J4TQYlUun0gpV7tgXoVC0YG4YjdhLzDIBbIoFIpziEsciIrzF6PRSFZWFlVVVedaFEUr0Ov1xMXFodPpHD5HKQNFs2RlZeHv70/37t2xxpcp3B0pJfn5+WRlZdGjRw+Hz1O5CYpmqaqqIjQ0VCmCToQQgtDQ0FZbc0oZKFpEKYLOR1u+M6UMFG5PVlYWV199Nb169SI+Pp45c+ZgMBgAePfdd5k9e/Y5lrApfn5+No9rtVqSk5NJSEhg4MCBLFy4EIvF0uxcx48f5+OPP24PMRuglIHC5Zw+fZrLLruMM2fOOD2XlJJp06bxu9/9jl9//ZX09HTKysp4/PHHXSCpbUwmU7vN7e3tTVpaGgcOHGDt2rV8++23PP1084m+HaUMkFJ2+M/gwYOlonNw8ODBVp9z7733So1GI++9916nr79u3To5atSoBseKi4tlSEiILC8vl0uXLpVTp06VkyZNkr1795YLFiyQUkpZVlYmJ0+eLJOSkmRCQoJcvny5lFLKnTt3ytGjR8uUlBQ5ceJEeerUKSmllJdddpmcN2+eHD16tFywYIHs1q2bNJvNUkopy8vLZVxcnDQYDPLIkSNy0qRJMiUlRY4cOVIeOnRISillRkaGHDZsmBwyZIh84oknpK+vr82/p/Hxo0ePypCQEGmxWOSxY8fkyJEj5aBBg+SgQYPkli1bpJRSXnLJJTIgIEAOHDhQLly40O64xtj67oCd0s59qZSBollaowz0er3EmrHa4Eev17f5+osWLZIPPvhgk+PJyclyz549cunSpTIqKkrm5eXJiooKmZCQIHfs2CFXrlwp77jjjrrxRUVF0mAwyOHDh8ucnBwppZTLly+Xt912m5TSqgzqK6+pU6fK//3vf3Xj/vCHP0gppRw3bpxMT0+XUkq5fft2OXbsWCmllFOmTJHvvfeelFLKxYsXO6wMpJQyKChInjlzRpaXl8vKykoppZTp6emy9j5Zv369vOqqq+rG2xvXmNYqA7VMULiMjIwMZs6ciY+PDwA+Pj7MmjWLY8eOtXlOKaVNZ1j94xMmTCA0NBRvb2+mTZvG5s2bSUxMZN26dTz66KNs2rSJwMBADh8+zP79+5kwYQLJyck899xzZGVl1c15ww03NPh9xYoVACxfvpwbbriBsrIytm7dynXXXUdycjJ33303p0+fBmDLli3MmDEDgJtvvrnVfyNYYzruvPNOEhMTue666zh48KDN8Y6Oay0qzkDhMqKjowkICKCqqgq9Xk9VVRUBAQFERUW1ec6EhAQ+++yzBsdKSko4efIk8fHx7Nq1q4myEELQu3dvdu3axbfffsu8efOYOHEi11xzDQkJCWzbts3mtXx9fet+nzp1KvPmzaOgoIBdu3Yxbtw4ysvLCQoKIi0tzeb5bfHgZ2RkoNVqiYiI4OmnnyYyMpI9e/ZgsVjQ6/U2z/nb3/7m0LjWoiwDhUs5e/Ys99xzD9u3b+eee+5x2ol4+eWXU1FRwfvvvw+A2WzmkUce4fe//32dBbJ27VoKCgqorKxk1apVjBgxglOnTuHj48NNN93E3Llz+fnnn+nTpw+5ubl1ysBoNHLgwAGb1/Xz82Po0KHMmTOH1NRUtFotAQEB9OjRg08//RSwPtH37NkDwIgRI1i+fDkAH330kUN/W25uLvfccw+zZ89GCEFxcTHR0dFoNBo++OADzGYzAP7+/pSWltadZ2+c09hbP7Tnj/IZdB7a4kB0NZmZmTI1NVVedNFFsmfPnnL27NmyqqpKSinl0qVL5XXXXScnT57cwIH43XffycTERDlw4EA5ZMgQuWPHDimllLt375ajRo2SSUlJsn///nLJkiVSSqvPoHZMLZ9++qkE5IYNG+qOZWRkyEmTJsmkpCTZr18/+fTTT9cdr3UgvvDCC3Z9BhqNRg4cOFD2799fJiUlyVdeeaXOUZmeni4TExPlJZdcIh977LG6OQwGgxw3bpxMSkqSCxcutDuuMa31Gbi8noEjqHoGnYdDhw7Rr1+/cy2Gog3Y+u7atZ6BQqE4P1DKQKFQAEoZKBSKGpQyUCgUgFIGCoWiBlcURO0ihFgvhDhU015tjisEUygUHYsrLAMT8IiUsh8wDLhPCNHfBfMqFIA1sq9+iK/JZCI8PJzU1NRzKNX5h9PKQEp5Wkr5c83vpcAhINbZeRWKWnx9fdm/fz+VlZWANeIwNlb9F3M1LvUZCCG6Y62UrNqrKVzKlVdeyTfffAPAsmXL6pKCAMrLy7n99tu5+OKLGTRoEF9++SVgrQMwatQoUlJSSElJYevWrQBs2LCBMWPGMH36dPr27cusWbM4F8F37obLEpWEEH7AZ8CDUsoSG+/XtVfr2rWrqy6r6EAefBDs5Oi0meRk+PvfWx5344038swzz5CamsrevXu5/fbb2bRpEwDPP/8848aN45133qGoqIihQ4cyfvx4IiIiWLt2LXq9nl9//ZUZM2ZQG/m6e/duDhw4QExMDCNGjGDLli2MHDnStX9cJ8MlyqCmFftnwEdSys9tjZFSLgGWgDUc2RXXVVw4JCUlcfz4cZYtW8bkyZMbvLdmzRq++uorXn3V2tGvqqqKzMxMYmJimD17NmlpaWi1WtLT0+vOGTp0KHFxcQAkJydz/PhxpQycnUBY8zbfBg5JKRc6L5LCXXHkCd6eTJ06lblz57Jhwwby8/Prjksp+eyzz+jTp0+D8QsWLLCb6uvl5VX3u1arbddSZ50FV/gMRgA3A+OEEGk1P5NbOkmhaC233347Tz31FImJiQ2OT5o0iX/84x916/7du3cD7Zjqe57iit2EzVJKIaVMklIm1/x86wrhFIr6xMXFMWdO0zCWJ598EqPRSFJSEgMGDODJJ58E4I9//CPvvfcew4YNIz09vUHxEkVTVAqzollUCnPnRaUwKxSKNqGUgUKhAJQyUCgUNShloFAoAKUMFApFDUoZKBQKQCmDDsFskVQZzZRWGSmvNmE0N991V9GUM2fOcOONNxIfH0///v2ZPHlyg/BiR9m0aRMJCQkkJyeTnZ3N9OnTbY4bM2YMF9r2t+qo5CLKq02cLq4kv8xAUaWR4gojxZVGKo1mzJamsRwaIdB5CPy9PAjy8STYx5NgXx0R/nrC/Dzb1J2nI/jb2tbfgM3x0ITeLY6RUnLNNddw66231jUqSUtL4+zZs/Tu3fL59fnoo4+YO3cut912GwArV65svdDnKUoZtJEKg4mM3HKyCis5VVRJcaWxVedbpKTaKKk2GsgrMzR4z9tTS0yQN3HB3vQI9SXY19OVonc61q9fj06n45577qk7lpycjJSSP/3pT/znP/9BCMETTzzBDTfcwIYNG1iwYAFhYWHs37+fwYMH8+GHH/L222/zySef8N///pd169bx/PPPk5qaWlcr4bbbbuPgwYP069evrnYCWBOh5s+fT3V1NfHx8SxduhQ/Pz+6d+/Orbfeytdff43RaOTTTz+lb9++lJWVcf/997Nz506EEMyfP59rr73W7jzuglIGraCkysiRnDKO5JRxuqgKSztFb1YazBzNKeNoThk/kEtEgBe9I/3pHelPoLeuXa7pztTe0I35/PPPSUtLY8+ePeTl5XHxxRczevRowHaK8h133MHmzZtJTU1l+vTpHD9+vG6uf/7zn/j4+LB371727t1LSkoKAHl5eTz33HOsW7cOX19fXnrpJRYuXMhTTz0FQFhYGD///DNvvPEGr776Km+99RbPPvssgYGB7Nu3D4DCwsIW53EHlDJoAYtFkpFXxr7sYk7kV3AuamDklFSTU1LN5l/ziAv2ZlDXIOLD/dx2KdFRbN68mRkzZqDVaomMjOSyyy5jx44dBAQEtDpFeePGjTzwwAOANV06KSkJgO3bt3Pw4EFGjBgBgMFgYPjw4XXnTZs2DYDBgwfz+efW7P1169bVLWcAgoODWb16dbPzuANKGdih0mBm98lCDmSXUFbtPumtWYWVZBVWYikv5KMXHuaLzz6hW9z5XQIsISHB5tq+ubyatqQo22v9PmHCBJYtW9bsdepfQ9poI9/SPO6A2k1oREmVkfWHc3h7cwY/ZhS4lSKoz+dvLyJtxzZuvu/P7DhecF7vUIwbN47q6mr+/e9/1x3bsWMHwcHBrFixArPZTG5uLhs3bmTo0KFtusbo0aPruifv37+fvXv3AjBs2DC2bNnCkSNHAKioqGhxF2PixIksXry47nVhYWGb5ulolDKooaTKyJoDZ3h3y3HSMoswmps+dUryc1j8yE2UFOQ6dLw9+HNqEg9P7MPW1cuQUrLpq48Z2iMUXx8ffs4sxHQeKgUhBF988QVr164lPj6ehIQEFixYwMyZM0lKSmLgwIGMGzeOl19+maioqDZd495776WsrIykpCRefvnlOqUSHh7Ou+++y4wZM0hKSmLYsGH88ssvzc71xBNPUFhYyIABAxg4cCDr169v0zwdzQWfwlxpMPPT8QL2nizCZGMLsD4rX1vAtm+WM/yqG5n+wAK7x0vyc3j/Lw9zy+N/IyAk3KXyluTn8NWSl9i3dR3G6ip0XnoSR0xg6l2PEhASTpCPjjF9IugR5prcfZXC3HlpbQrzBeszMJkt/JxZxM4TBVQbm3+a/jk1CZOhuu711tXL2Lq66dqv9rjQaEBK1nz4egOl4QoCQiPQ+/hhMlTj4emFyVCN3sevTukUVRhZtTubnuG+jOkdQaDPhbf7oGgbF+Qy4WhuGe9vO8GWI3ktKgKAJ95bR8rYVHRe1hp6Oi89KeOm8Mg/VzU4Xou0WJBSsnX1Mh6e2Ic/pya5VP7SonwuTZ3BnEWfcGnqDEoL85qMycgt5/1tx/npWAGWFiwehQIuMMugoNzAD+k5HM+rcGh8fXPf1tM4Nr5fk+MAGo0Wi8XcwIRvzbVaWlrcNv8359S198+3O85kkWw5kseRnDImJkQS5udld6xC4RLLQAjxjhAiRwix3xXzuRqT2cLWI3l8uP2Ew4oAYM1Hb3Bs/07WfPi63adxaVE+CNFgGWGxWAtvGqurGpjwtdhyONZea/Vbr7L4kZvIPnrIZU7JsyVVLPsxkx3H22YlqAYjnY+2fGcucSAKIUYDZcD7UsoBLY3vSAfiyYIKvj90lsIKx8OFG/sIavHw9OLl1XubHG/s1BMaLX0Gj8AvKJTqirIGT3Jo6HD8ac3nNq8FVi96Y2dlW6m1PP780hvcOGagw5GMx44dw9/fn9DQ0As+yKmzIKUkPz+f0tJSevTo0eC9dncgSik31rRWcxsqDWZ+SM/l0OkmzZ1a5In31jW4uT08vdB7+3L3S0sbjLO3jDAbDYRExja5iW05IgGE0CBlU99Frd9h6+pldhVRY+wtN2otj/ffWIjZ+xkm9IukV6R/i/PFxcWRlZVFbm77b5sqXIder6+LwHSUDvMZdGR7tV/OlPDD4VwqDG2rk2/LY19mqGbb6uUNbvDGy4hhV95C78G38uN3u8nY78fK1yKoKNVSXSmortQQHptHcV42lWXZSFmM0FQQ0cULve8pThxcARwBchrI4uHpRdLIiYy97g8sfuSmJjd545u//nKj4Gw2mb/swWT8LRGqvnL5Lu0Eo3uHo9PaXy3qdLomTxfF+YnL4gxqLIPV53KZUFJlZP0vOWTkljs919KnZ7N/2/dIS3O7DWHAyHo/KcBv5rePvxm/IBNe3hJPvQVPvYXTx45SlFsGIhCkH0LEIWV9nZwNbAR+ALER5C9cmnojQLMxDghhW1YhSBlzld24hDB/L6YkRRPkc2FnRl4oNLdMOC+UgZSSPVnFbDmSh8Hkugg8WwE+vVNmER47hx1rtZQX98Xqg63CL/AoyWMC6d5fS1iMgbAYIybDmSYm+9KnZxMQEs6wyTew/dsVFOfnM/Wuf5KbpWP1219jrB5IYU4vzKbap38m8AmwAnDuM6tdwjRWKF46DZMSoogPd590WkX7cF4HHeWVVfP9obOcKqpq8xy1pvY1f3ycL954vu7mrV0uGKu90GjvwFh9Kwe2WT9HH/8TwLNoPdZjNm0nadQ0pt23oMG8K1/7bRlRe/PZ3hY0EhZjpN/QywGQspC8U2Uc2eND2g9+HEmbg5RzgaNotMvpP+xXJt10Des/eauBogoMjSTv1IkGPgidlx4f/yB6DRrOZdN+z/ZvVzTZoag2Wvh6zymGdg9heLxyFF6ouGo3YRkwBqvdfBaYL6V82954V1gGJrOFn44XsPN4oc1KQq2h1tSO6BpPTuZRBl9+NQVns7n8xn+zctEpivPHYjHr8Q04jn/IGm5fMIGv/313gyd8SUFu3Y3e2t2Illj+6sv8tEYixEykHIsQFgaNqcBoeJn9W/+CVueJ2WggODKWfhePprSogL2bvkMIDSBbtSPRI8yXKwZEoddpWy2nwv3pkGVCa3BWGRzLK2f9Lzmtri7UGHs3LUwC5gGX4am3MGhsKcMnF9GldzVCtBwg1FL+QGupv7RY8+F3HNqRiEZ7J4ZKD/yD07niljxOZbxep5AaL0XqKypHCPH1ZOrAmAu+wtL5yHnTXq240shXe06xane2U4qgNuhnzqIVpIxNxcPTCxDAdGAX8B3QE3gQQ1UQu76PpGsfqyKAhrsItrCXP4CUDQKJHM12vG3+Yq69fz6x8X3xDy7CbLyPQZfdxTV/zEGr7cGniy4l5+S/GXvdW03GX3v//FYpArBGai7bkcmJfOcdsYrOQ6ewDKpNZnYeL2R3ZqHN1OLWUj/oR1ok274tBv4CJAOHgReBj9B5aRs80Vtj/tt6OvsHhzXYEbCXBWkLe9fW6gKYcsevrFsWQlmRBwnDyph6dy7hsc5ZTWAt2jqqdxgpXYOdnkvhHnTaZYLZItmTVcRPxwqobGPMQH2a3lBDsd74Y/HwzMI3YCEWyzJKC87YXG+31fz/81WJDfb67WFLqdR3bjZ2GNa/dnWlYNOqIL5fEYLJKBg7vZDxMwrw1Dv//SbFBTK2TwQajXIsuiNHckq5KKLlADLoxMuEZT9l8sPhXJcoAvgt+9DDsyfwIfAjHrpBXPn7o/zliwrmf3wP3fslMmLKTB5+44smGYEtpQ/bI/myyYA1gQmsN31wZGzN8uS3LMgn3v++ybm1S5Jt36xo9tpe3pLxMwqZ9/ZxkkeXsW5ZKC/d0Z29m53fLtybVcyqtGyqjK75HhSuY392Md/sPeOS/BG33lqsMDhecsyRrD8f/0jyTs3EZJiC9U9/nsGXH2XCzMfqxrSUEVibsFTf/LdHY0ukNoHJZKjG08sbs9FgV6nYDV3WaJiz6BO71w4INTPr0TMMn1zMZ4sjePeZGAaOLuXa2Tn4BbX9Zj6RX8EnO09y9cBYVSPBTdhxvIDNvzZNX28rbr1MWLLxKOXVjv0Hbm79XZKfw5vzllNd8SKFOX4Ehf/ItfeX8suOpa32tLeG5hKY0ndtJnHEBLsef1fsSJjNsP6TEP77QSjevmaufSCHgaPKnPqbfDy1TE2OITrQ26l5FG1HSsmmX/PYdaKw7tiD43s5FB9yXgcd2atCVLv+Li3Usngu5GUvRutxDJhO/0tCSBi2gC697uX9vzxMSUGuy8uTQdNlhb0EJlsWSGuXJLYsI60Wxs8oIGFYGctejeK9Z2MYNKaEa+/Pwce/bZGaFQYzK3dmccWAKIcSnRSuxWKRrDl4tk0JeC3h1j4DR7BXhejx975n7pWvMv+GIPKyhwLPYDb1A/5bV4Ho6VmXNbtF6AocqUrkinOb2+6M7mFgzqJMrrw1jz2b/Pnrvd04dkBvYxbHMFkk3+w7zc7jBW2eQ9F6DCYLX+7JbhdFAOfJMmHlovls+3ZFXSTeoLH3U5jzNMf2B+ETsI/qypsxG/e0OE9bIwTPJS1tdza2GE4c0vPBC1EU5eiYdEs+l99QgMaJYMPE2EDG9VU7De1NhcHEl2mnOFNsO+zeFcuETmcZ2ArUqX2CPvC3T+jS+1/8/L/nOXHIE7gXv8BrsZj2WouUQt2/YTHdmlgTtrz57o49y6j2b2lsMXTrV8Uj/8xk4OhS/vNuGK/PjeDvDzzQ5opK+7KtOw3VJrXT0F4UVxj5ZMdJu4rAVXQ6n0H9/9z1k39ys3WsWBhJ5mEf4L9YzHcCJ8k5+du5Ud164RccSmSXnhza8UOrtwibQwhrGG+Ir7WjcpCPjiAfT7x1WnRagU6rwVOrQWI19wxmCwaThbJqE4UVBooqDBSWG8kprW7VFp4938Jzt1zerC/lpnln6J1Swad/D8ZieZuVry3h9gWzHLqmrcSu8moTl0RpuPP3N7NixYo29y9QNORMcRVfpmW3uTZHa+g0ysCeo1Cr0zP1rhN88bovYABmAw0rEqWMm9LEC7/06dn0u/gyh7YI7RER4EWXYB9igryJDfLG29Mxe9vbU4s31rHh/l704LceB1JK8soMZBVWkF1USWZBRYsVnG1tdzau1tS4OOujqbWBUEnA5+zf+jgPT3wIjcdbdO83sNkt2lqF/OGLc8nJPFqnmO986Gk2bd7MM888wxtvvOHQZ6Gwz9HcMv6z77RLom4dodP4DOzVFigveZXjB4KITyrE2/chDv+8AmO11Zyq7V/gqjqCQkB0oJ6LIvy5KMKvQzoim8wWThRUkH6mlIy88lbVa2jsS6n/OXz88qPsXLeqppKzP0LzIdJyFSHRGyg4PZlLU3/XYtm2ltDr9Q1amyscZ8/JIjYcznW40/cFsbVoq86gVueFsfpmDv20CJ2XJ9c/dIZLrijhs9eqMRmq6/L5k0ZOxC8wxOkKw35eHiTEBjAgNpAAfccG3HhoNcSH+xEf7ofRbOHXs2WknSzibEnL60dbFkPj0GhrIFQR0jIFeJyC008DG9m6+ndsXd0Hrc6TV76xthavtTb2blnbQClodZ4EhIRTWpiHyVCNzkvP5VeksvTNf7j40zj/sRVD0FG4vTJoXGdw8OX3kXfqSY4fjAB+4N6XfOnaJwCw/Z+/ub4CLdE1xIeBXYLoGebrFt5ynVZD/5gA+scEkF1USVpmEUdyyuw+PWxFUza0CMx1gVA5WccoOP0cQuxFyg+AHcA1DBoTXTdHrX/CbDQgNBqkxYLQaLCYjE0iKkvMOvbkQ3iEROsGn11nwGCy8N2BMxzNaV1gWEl+DmPG3Om0r8Ztlwne3t5UVTV++t0CLMLTK4DYXu9xbP8dXJp6g0tbmAkBF0X4MbR7CBEBbd+L7ygKyw1sz8jn8NlSmvsqW2fi9we+AuKw1rB9v87xWJuNefZkBmWF+XUO2X1b1tqMqIwN9iY1KRofT7d/7pxTSqus6fk5Jfa/I3sh9ytfW8D2b1dw9913t+ir6ZRZi2lpaYwaezlVFeWYDCHAv4ApwCbgNuBog/H1zdm2IAT0jQpgaI8QQjphUY+8smq2Hs23+1SxFxrt6e3LiYO7qSgtsvpahAb/4FBKCwG5DLicyK6fcfeLUQSFtW23xV/vwdTkGCL83V+5ngvOllTxVdopyqqbz8WpDbmvrcTVuPJ1Lc35ajplnMGSJUsoKyrAZLgOOACMBx7GWl3tN0UgaiJmBo2Z3OZr9Qz35aZh3bhiQFSnVAQAYX5eTB0Yw/UXdyHcv2kbtcZbkEgLIZGx/P6Jv5NwyZi64wKJj18ggjy0uquBNzmbeS1fvJ5EXnZem7o8lVaZ+GTHSQ6fKXXRX3v+cOh0CZ/uPNmsIvhzahIPT+zD1tXLkFKyc90qMvbtwGQykjI2Fa1uJDAOHx8fZs2axbFjx9oki6vaq10hhDgshDgihHis5TPs4+3tjRCCf/5zFVZT9QOsymAg8Dd0Xp6ExXSrGy9rMgF3rl3V6iansUHeXH9xF65Ojj1v+hDGBnkzc2hXxvaNwEvX8OttrkVc/eOVZcVcmjqDB1/7mCHj16Dzmse+rX4sejCKjH3ZbQrfNpol3+47zeZf81S7Nqw5Bj+k5/Ld/jMtbh3WBpY1QUp+Xh+B2bgOIV6hsrKKgICANvsNnF4mCCG0QDowAcjC6nmaIaU8aO+c5pYJp0+fZu7cuSxffisWyyg02vkER6yk4MzxBoU/4xMvprQon/Sft2Axm1uV1eev92BkrzD6RgW0+e/uDFQYTGxMz3Mqln3lawtq0qevApYDRUAqsKfN4dtdQ3y4MjHqgvUjVBnNfLP3NJkFjvf9XLloPlu/WV6v8rUOjWYxFstd+Ifs4avPNHz66ZucPn2azz//3O487b21OBQ4IqXMqLnYcuBqwK4yaI7o6GgCAgKQ8n48dB6YTYeQMrbJLsGMP73IykXzkRaLw1GEWo0gpWswQ3uE4Onhtiskl+Hj6cEVA6LoHenH94dyWlyT1qepw/EbYASwGthEfNJCbv6/8W2SK7Oggo9/zOTKxGhigy6sVOjTxZV8s/c0pVWOfxdgtd5GTJlZU/l6D/ApFssIPPWv8dDiwVx22aWMGeNcwp0rLIPpwBVSyjtqXt8MXCKlnN1oXP32aoNPnDhhd85p06YRHR1NyOArWb9qmd2aA62pAtwlxIfL+0ZcsBV/q4xmNhx2vPekLYejtJjRenTHbPoMSOba+3MZMaW4zTJphGBkr1BSugZfEL0adp0oZMuRPKdK+78+9zVOpj+DlAH4+D9Mcd4/uDR1Bpu/+sjpoCNXKIPrgEmNlMFQKeX99s5pj+Im9tDrtIzqFcaA2ECn5jlfOJJTytqDOQ7lP9SPYDQZqonsdhE3PfZXtny1igPb76K0cBhjpheQekceGicMre5hPkzsH4Wv1/m5bKgymll78CxHWhk/0Jif1/uz/K+RmAyZWHfW9jd435GIz/beTcgCutR7HQeccsG8TtMnyp9bL+2mFEE9LorwZ9awrsQGt2ye13csjpgyk4i4HsTG9+X6hx7jqY9DGDG1iA0rQ3jv2WgMVa1/stdmoO5NP8EH2084fbO4IycLKviwlX9b48zcotwcnpm1iQ9fiKZr7yr+9K8zpIzt3iBTdebMmW3eRajFFcpgB9BLCNFDCOEJ3Ih1G+Cc4eOpZcrAaCYnqmAXWwTodUxPiWNYz1A0zZiWzfVf0Gph2n05XH1PDvu3+vH6n7pQUtC6wgj1o0srDWa+3nOKdQfPnhfp0EazhfWHc3hn7c+8cN+NNrdj7fXNqP+5VFcKFj/iSVHuH4josp57XsoiukdIk0xVZ3YRanFVe7XJwN8BLfCOlPL55sa35zKhd6Q/4/pGOJxBeKFzsqCCb/eddipFdv9WXz58IRrfIDN3PptNVPfmy8K3VJDFX+/BmD4RXBTRORvBni6uZM2BsxSUG5qtzdlyEFEM8DXWbfW5WG8x6+fU7+LRDfxlER6Vze4i1NIpIxChdcpAr9Myrm8EfaJUXb7WUlpl5Ju9pzntRPGMk+levPVULMYqwa1PnqbPYPvbZo4We70owo8xfcLx7+DksLZSZTSz5Uge+7KL+dNV9hUeYDs0XAhSxlxV87kkAF+B8EejvRmL6ctmt88vyEpHtugS4sNNw7oqRdBG/PU6rhvSheQuQW2eo0vvah58LZPgSBP/fjyWrat/89M0NocdLfZ6JKeM97edYMfxAkzmthVw7QiklOzPLubdrcfZm1WMlM1XoGo+iGg1xuqrgI2AgZDIG5Hmr1xWhKc5OvWCWqsRXBofyuBuF8bWVHui1QjG9o0gMkDP94fOYmrD9ldwhIn7/5bJhy9Es/K1SHJOejL1rty6NfDqt16l4Gw21/zxcdI2fcfg8b+z2ya+FoPJwuZf89hzsohhPUNJiAlwq+/6RH45W4/mNylJ1pLC0/tYl0C/BRGBh6cercfTVFf8meieBcRd9BKHd+5wuE+Hs3RaZRDko2NyYjSRnSCzsDPRPyaAYF8dq/ecdjhIqXE23fQH0lj0QCYbv5jJxi+2YA1UssbUA/z13t8hhMDTS1/nnGyJ0ioTaw+eZXdmIUN7hNIrwu+cppWfLKhgW0Y+2YX2t/Kaa7jTMIjoO8AHk+FtTIYbGTK+hOsezKeydBb5p3cxYdYfCQgJdyod3xE6pc+gd6Q/4/tH4OWhnITtRVm1ia/32K/GW5/GTrLa1136vEHm4TtB/gJMBTJsnt+WsGZ/vQfJXYIYEBuIXtcx/w/MFsnR3DL2nCwiqxkl0BqWPj0bL59enEx/jrMnAonq9g5/WjIaIZpvDNQYV/gMOpUy8NAIRvcOZ6ATa1uF45jMFtYePMsvdrINHauRMBb4FGvL+xuAdXXveHh6kTRyYqu6RDXG00NDrwg/ekf60yXEp10KqRRVGNifXcKBU8UuL0x64pCepc/EUFWu4aZ5pxkwvLxV3b5ruaAciIHeOm64uItSBB2Ih1bDlYnRDI8Ptfl+YyeZrYaygWGHSBzxBHqfYuA74EHAWp/SbDQ45BCztx8PVp/CgVMlfLE7myUbM1hz4Ay/ni2ltKrtLenNFklmfgU/pOfy/rbjLN1ynB3HC1yuCH78LoDFc+Pw0Fl44O+ZDBheDrRc/r696BQ+g+5hPlw5ILrDzEFFQ4b1DCXYx5M1B840cCzaah/XuPxZwrCxTH/gId56ag5njj9JwZm/ofcdQ0zPvxPdPcYhh5it8vi2qDKaOXCqhAOnrPkX/noPIgP0hPl54eulxcdTi4+nB14eGswWidEiMZsl1SYzBeUGCisM5JcbKCw3tGtFYpMRvnwzgi1fB9E7pZyb/+80vgG/7Za0tdu3s7i1MhAILukRwvD4ULfyIF+I9Inyx1/vwdd7TjV4QjZ2ku3bstam0+yOZ17FYoHvl+Xx3ftTqSi5gpFXnyKya9MneK1DsnEQTuPeDy1RWmWitKrMrcKci/O1vP98NMf2+zBmegFX/SEPrY1nXGu6fYf6ebrk/nBrn8HJggq6hPh0gEQKRymuNPJlWjb5Zc1HGTZH+s8+fPBCFCaDhusfOsOgMQ1v1vqReRazyalO1O7E4V0+fPRiFIYqDdc/fJaUsQ19MfZqHNpDCBjcLZjhPUPx0Dq24u+0PgOlCNyPQG8d1w/p4tR30zulgkfeyCS6RzUf/CWGla9FYKgSNst7WYNwqjrUXHYGW/4NixlWvannX/Ni8Pav5qHXM5soAmi+eW5jgn2s38OoXuEOK4KWcGtloHBP9Dot1wyKpX9M2ytFBYWbuO/Vk4yZXsDW1UEsvK8rNz78H3wDQxo5IKMYMuGaNnWxPhc0vqGLcj1487E4Nn7eFXiPngNm4+2b1UBhNFaCtV3C7ZXwS4oLZNawbsS4uDCMWy8TFO7Pjxn5bD2a79Qc6bu9WfZKFCX5AimfAF7Bw1PXpAuUO2N7O3Am8Dqgw9r2790G716aOoPpDyxwOFfD00PD5f0inCrX12mXCQr355KeoVyZGIWHE/v7bz3Zi+K8SKT8AngR2ITJEA9CuL0lULssmLNoRb3twFCEWAl8RGS3avpf8hA6r+VNzq21AJ67dXyLuwdhfp7MGNq1Xet2KmWgcJq+UQFckxLb5rRx6776CDw8bwFuBnojRBqjrj7OTfOcq+vX3tQuC7Z9swIvbz+M1dOA/Ug5BXiUngMeJiistKYtoLXknqamvH/9+AF7lasB+kX7c+PQru1exl8tExQuo7DcwKq0bIoqWh/w07DEWgDhsZ+Rmz2K8FgD02bntJgS3dgL31rPvCNz1qfpsqAf1iXBWKz1fu4Aftv+FBoND7/+BR+9+AhnThypi8tobhmkEYJRvcNI6RrcavntoZYJig4h2NdqysbVlFRrLnKwMQ1LrE0gqvuz3Px/+yjMPcu/5sXxr3mxZB+13duivtOu9pqr3/mrTc+8ozK15NmvjRL08IwCXgL2oPUYwviZaQwa8xw6r3Tgt6f//I83Ehvfl/C4HoyYMrNFh6i3p5ZpKbEuVQQtoSwDhcsxWyTrDp3lmccedjjRxhbWng2f073/W+ScnEFlmYaUcaVccUs+odHGVvWPrA1Uain5x96cGg8d3fsl11kKOSfzWTh7I4bKu4Eg4G0unvATM/70SAMrpy1O0HB/L6YMjCHQ2/VFXTptopKic2K7aa7j2Ym2b8hAhOZxtB6PYDYJkkaUMWR8Brs3PFHnhbdHyrgp7N30X5t9CWtv8mv++DhfvPE8E2bdy4cvzKW6shyTobquRHx4XHfysk8wZMIdhMX8hbUf6TEZ/fDyXsP1D3uSsfffdaX6W1PCvzG9Iv2YlBCFzkWxA41pN2VQUyZ9AdYF01AppUN3uFIG5ze1XbFWrVpFRUVFqyMH7W21jb3uD3y6aAldev2Ln9dHU1mmRaPdicW8CI3Ht1hMBXWt5sG6TkdKhl91IxNn/dFmH4jamzyiazw5mUeJ6BrP2RNHbEg1Gmvbj+mAF9aavwuA3YC18W+3vgPb7KMQwpoDckmPkHYNvW9Pn8F+YBrWGk0KBfBbV6yqqir0ej1mo4GAgAC7N4mjZdG2fbOCk4f/g5R/5qmPMuiR8C4WcyDwAYIcvLz/i8VyPWAtuZY0cmLdurx2zloLorZHZ27WcaSUnD1xpO5fKzqsCuAF4BfgB+AqhOYtAkLG4+F5PbC7zicwaMxkh6MHG+PpoSE1KYZhPc9tDo6rqiNvAOYqy0BRS21XrLvuuoslS5aQlX2Km554zWZREFvr+Pqm9sL7rkFa7NVAFMClwHVYn9qxCCHx9juJl89+JsxMITa+moBQE5/94z68fX0oLcrn8K7N9eb0BHphNXD7AoOBcUAAYAQ2A++C+AxBBRFd4sk5ebSuuYwtHF0SBfvomDIwhtAOavzb7j4DR5RBa9qrKc5PajsPp50sAloumV5L42WDh6cX/sFhlBbmYTJU1y0jUu94lMKzXTi8y5cTv+jJPKynsvS32AchJL4BZkzGXKoqjIAe8AZ8+c1ItmCtyPQ9fkE7SBgG+WcOUFaYj19wKJFderJvy1oSR0xg2OQb2Pj5u6Tv3kZFaVGrk6l6hvsyKSGqQ1PznWq8KoRYB9jqzvC4lPJLR4WQUi4BloDVMnD0PMX5g6am6GpUoLXo6hPvrbMbhlsfR+om6H38CAoLJyisih4J1qWAxQIfvfgvdm/IoGdiKr2Sr6a00IMD2w8REKqhqiIfk6EEoyEHH/8zVFftIjC4CP9QXyK79KSkIJcbHm7q+Ktfi7C2AfC2b1c4nEwlBAztEcLwc7wsaEyLykBK2bZWuwqFHfpFBxDu78U3e3UOF/Gon9+/8fN32bX+62arKze2OjL2rSJj3x2NrI4Yl/w9rak9oNdpmZgQSXy4+zWIceviJorzlzA/L2YM7cq/qoodupHqb83pvPRIs7nZ6sqOWh2uoL5szVUwjg7Uc2VidLvED7gCp5SBEOIa4B9AOPCNECJNSjnJJZIpzns8PTRsXfct+7OL+SE9l9j45kuBN37aN1f56FyVDrPHoK5BjOoV3i4FW12FU1uLUsovpJRxUkovKWWkUgSKtjAgNpBZl3QlOrD5HhitLRTaXPJPW2lNiDWAr5eWqckxjOkT4daKANQyQeEmBPl4cv2QLmw/ls+OY4VYbOxytfZp76j53hocLc4K1l6R4/tFdpomwEoZKNwGjUZwaXwYF0X48f2hHJsNXFrjrHMlrVmieOk0jOkd4VQlqHOByk1QuCVSStJOFrH1aD4Gk/NNV12R0uxINaK+Uf6M6h2On5d7PmdVCrOi0yGEYFDXYG4e3o2+Uf44ux3fmmKjtmhpiRLq58n0wXFcmRjttoqgJTqn1IoLhgC9jisToxncLZhNv+aRWWC/yIktWmPet4StJYqPp5Yh3UNI7hLk9g7CllDLBEWn4kR+Odsz8jlV1HJDWHDcvG8tPp5aBncLZmCXoHZLN24PnApHVijciW6hvnQL9eVMcRU/Zxby69kymzsPtbg63iDIR0dSXCCJsUF4enQeJeAIShkoOiVRgXomJ0ZT0svIgewSjuSUkmeny5OzOxBajeCiCD8SYwOJC/Z2q3wCV6KWCYrzhvyyatLPlnE8v5zc0mrMlrb/3/bz8qBbqE+NJeJz3jT9VcsExQVBqJ8Xw/28GB4fislsIae0mjMlVeSVVlNhMFNuMFFRbabaZEar0eChEWg1Ak8PDcE+ngT76gjx9STcz6vD6gu4E0oZKM5LPLQaYoK8Xd6C7Hzm/PKAKBSKNqOUgUKhAJQyUCgUNShloFAoAKUMFApFDUoZKBQKQCkDhUJRg1PKQAjxihDiFyHEXiHEF0KIIBfJpVAoOhhnLYO1wAApZRKQDsxzXiSFQnEucLYg6hoppanm5XYgznmRFArFucCVPoPbgf+4cD6FQtGBuKS9mhDiccAEfNTMPPV7LbZJWIVC0X443V5NCHErkApcLpvJh1a9FhUK98bZjkpXAI8Cl0kpW1ecTqFQuBXO+gwWA/7AWiFEmhDiTRfIpFAozgFOWQZSyotcJYhCoTi3qAhEhUIBKGWgUChqUMpAoVAAShkoFIoalDJQKBSAUgYKhaIGpQwUCgWglIFCoahBKQOFQgEoZaBQKGpQykChUABKGSgUihqUMlAoFIBSBgqFogalDBQKBaCUgUKhqEEpA4VCAShloFAoanC2vdqzNa3V0oQQa4QQMa4STKFQdCzOWgavSCmTpJTJwGrgKedFUigU5wJn26uV1HvpC6h+CApFJ8Wp6sgAQojngVuAYmCs0xIpFIpzgmimCZJ1gAPt1WrGzQP0Usr5duapa68G9AEOOyBfGJDnwLhzibvL6O7ygfvL6O7ygeMydpNShtt6o0Vl4ChCiG7AN1LKAS6Z0DrnTinlEFfN1x64u4zuLh+4v4zuLh+4RkZndxN61Xs5FfjFmfkUCsW5w1mfwYtCiD6ABTgB3OO8SAqF4lzgbHu1a10liB2WtPP8rsDdZXR3+cD9ZXR3+cAFMrrMZ6BQKDo3KhxZoVAAbqIMhBBXCCEOCyGOCCEes/G+EEK8VvP+XiFEipvJN6tGrr1CiK1CiIEdKZ8jMtYbd7EQwiyEmO5u8gkhxtSEth8QQvzQkfI5IqMQIlAI8bUQYk+NjLd1sHzvCCFyhBD77bzv3H0ipTynP4AWOAr0BDyBPUD/RmMmA/8BBDAM+NHN5LsUCK75/cqOlM9RGeuN+x/wLTDdneQDgoCDQNea1xHu9hkC/we8VPN7OFAAeHagjKOBFGC/nfeduk/cwTIYChyRUmZIKQ3AcuDqRmOuBt6XVrYDQUKIaHeRT0q5VUpZWPNyOxDXQbI5LGMN9wOfATkdKRyOyTcT+FxKmQkgpXRHGSXgL4QQgB9WZWDqKAGllBtrrmkPp+4Td1AGscDJeq+zao61dkx70dpr/wGrdu5IWpRRCBELXAO82YFy1eLIZ9gbCBZCbBBC7BJC3NJh0llxRMbFQD/gFLAPmCOltHSMeA7h1H3idG6CCxA2jjXe4nBkTHvh8LWFEGOxKoOR7SqRjUvbONZYxr8Dj0opzdYHW4fiiHwewGDgcsAb2CaE2C6lTG9v4WpwRMZJQBowDogH1gohNsmGCXvnEqfuE3dQBllAl3qv47Bq3taOaS8curYQIgl4C7hSSpnfQbLV4oiMQ4DlNYogDJgshDBJKVe5iXxZQJ6UshwoF0JsBAYCHaUMHJHxNuBFaV2gHxFCHAP6Aj91jIgt4tx90pFOGjtODw8gA+jBb46bhEZjrqKhY+QnN5OvK3AEuNRdP8NG49+lYx2IjnyG/YDva8b6APuBAW4m4z+BBTW/RwLZQFgHf9fdse9AdOo+OeeWgZTSJISYDfwXq0f3HSnlASHEPTXvv4nV+z0Z6w1XgVVDu5N8TwGhwBs1T16T7MDEFgdlPGc4Ip+U8pAQ4jtgL9bw9reklDa30M6VjMCzwLtCiH1Yb7hHpZQdls0ohFgGjAHChBBZwHxAV08+p+4TFYGoUCgA99hNUCgUboBSBgqFAlDKQKFQ1KCUgUKhAJQyUCgUNShloFAoAKUMFApFDUoZKBQKAP4fc1i+G5OUWJ8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "class SimpleSincKernel(gpytorch.kernels.Kernel):\n",
    "    has_lengthscale = True\n",
    "    \n",
    "    # this is the kernel function\n",
    "    def forward(self, x1, x2, **params):\n",
    "        # apply lengthscale\n",
    "        x1_ = x1.div(self.lengthscale)\n",
    "        x2_ = x2.div(self.lengthscale)\n",
    "        # calculate the distance between inputs\n",
    "        diff = self.covar_dist(x1_, x2_, **params)\n",
    "        # prevent divide by 0 errors\n",
    "        diff.where(diff == 0, torch.as_tensor(1e-20))\n",
    "        # return sinc(diff) = sin(diff) / diff\n",
    "        return torch.sin(diff).div(diff)\n",
    "\n",
    "# Use the simplest form of GP model, exact inference\n",
    "class SimpleSincGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super().__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = SimpleSincKernel()\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "# initialize the new model\n",
    "model = SimpleSincGPModel(train_x, train_y, likelihood)\n",
    "\n",
    "# set to training mode and train\n",
    "model.train()\n",
    "likelihood.train()\n",
    "train(model, likelihood)\n",
    "\n",
    "# Get into evaluation (predictive posterior) mode and predict\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "observed_pred = predict(model, likelihood)\n",
    "# plot results\n",
    "plot(observed_pred)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
